Last updated: 2024-04-11
Checks: 7 0
Knit directory: DGRP_sexual_conflict/
This reproducible R Markdown analysis was created with workflowr (version 1.7.1). The Checks tab describes the reproducibility checks that were applied when the results were created. The Past versions tab lists the development history.
Great! Since the R Markdown file has been committed to the Git repository, you know the exact version of the code that produced these results.
Great job! The global environment was empty. Objects defined in the global environment can affect the analysis in your R Markdown file in unknown ways. For reproduciblity it’s best to always run the code in an empty environment.
The command set.seed(20210706) was run prior to running
the code in the R Markdown file. Setting a seed ensures that any results
that rely on randomness, e.g. subsampling or permutations, are
reproducible.
Great job! Recording the operating system, R version, and package versions is critical for reproducibility.
Nice! There were no cached chunks for this analysis, so you can be confident that you successfully produced the results during this run.
Great job! Using relative paths to the files within your workflowr project makes it easier to run your code on other machines.
Great! You are using Git for version control. Tracking code development and connecting the code version to the results is critical for reproducibility.
The results in this page were generated with repository version 3181c8d. See the Past versions tab to see a history of the changes made to the R Markdown and HTML files.
Note that you need to be careful to ensure that all relevant files for
the analysis have been committed to Git prior to generating the results
(you can use wflow_publish or
wflow_git_commit). workflowr only checks the R Markdown
file, but you know if there are other scripts or data files that it
depends on. Below is the status of the Git repository when the results
were generated:
Ignored files:
Ignored: .Rapp.history
Ignored: .Rhistory
Ignored: .Rproj.user/
Untracked files:
Untracked: %
Untracked: Antagonism_trait_plot.pdf
Untracked: Chapter_2_Figure_S1.pdf
Untracked: Chapter_2_Figure_S2.pdf
Untracked: Conflict_plot_females.pdf
Untracked: Conflict_plot_male.pdf
Untracked: Diff_atlas_plot.pdf
Untracked: Figure_3.pdf
Untracked: Figure_4.pdf
Untracked: Figures.Rmd
Untracked: Mainz_response_space_plot.pdf
Untracked: Manuscript/
Untracked: Morrissey_model.html
Untracked: Morrissey_model.qmd
Untracked: Mutation_load_phenotype_effect.Rmd
Untracked: PRISMA.pptx
Untracked: R_atlas_plot.pdf
Untracked: R_multivariate_brms_simple.rds
Untracked: R_trans_multivariate_brms_simple.rds
Untracked: R_transcriptome_medians.csv
Untracked: Reported_heritability.xlsx
Untracked: Rplot.pdf
Untracked: Rplot01.tiff
Untracked: Useful_cuts.Rmd
Untracked: add_later.Rmd
Untracked: antagonism_R_space_plot.svg
Untracked: antagonism_trait_plot.svg
Untracked: big_f_mutation_plot.pdf
Untracked: big_m_mutation_plot.pdf
Untracked: code/get_gene_annotations.R
Untracked: conflict_plot.pdf
Untracked: data/DGRP_mutation_load.csv
Untracked: data/R_medians.csv
Untracked: data/X_burden.csv
Untracked: data/autosomal_burden.csv
Untracked: data/female_X_mutation_corrs.csv
Untracked: data/female_autosomal_mutation_corrs.csv
Untracked: data/female_mutation_corrs.csv
Untracked: data/input/
Untracked: data/male_X_mutation_corrs.csv
Untracked: data/male_autosomal_mutation_corrs.csv
Untracked: data/male_mutation_corrs.csv
Untracked: data/organismal_level_output/
Untracked: data/transcriptome_output/
Untracked: evidence_ratios.R
Untracked: female_mutation_load.pdf
Untracked: figures/
Untracked: fits/
Untracked: height_example.pdf
Untracked: mag_R_sem.pdf
Untracked: male_mutation_load.pdf
Untracked: mutation_load_plot.pdf
Untracked: sem_figures.R
Untracked: sim_me.rds
Untracked: sim_no_me.rds
Unstaged changes:
Modified: _workflowr.yml
Modified: analysis/Transcriptome_analysis.Rmd
Deleted: code/export_gwas_results_for_shiny_app.R
Deleted: data/all.dgrp.phenos_unscaled.csv
Note that any generated files, e.g. HTML, png, CSS, etc., are not included in this status report because it is ok for generated content to have uncommitted changes.
These are the previous versions of the repository in which changes were
made to the R Markdown (analysis/Main_analysis.Rmd) and
HTML (docs/Main_analysis.html) files. If you’ve configured
a remote Git repository (see ?wflow_git_remote), click on
the hyperlinks in the table below to view the files as they were in that
past version.
| File | Version | Author | Date | Message |
|---|---|---|---|---|
| Rmd | 3181c8d | tomkeaney | 2024-04-11 | Completed analysis |
| html | 24fdbb4 | tomkeaney | 2024-04-11 | Build site. |
| Rmd | 1228613 | tomkeaney | 2024-04-11 | Completed analysis |
| html | 1228613 | tomkeaney | 2024-04-11 | Completed analysis |
| html | ad196cf | tomkeaney | 2024-04-11 | Build site. |
| Rmd | ec2fc6d | tomkeaney | 2024-04-11 | Completed analysis?? |
| html | 5b2cdee | tomkeaney | 2023-10-19 | Build site. |
| Rmd | 4188afb | tomkeaney | 2023-10-19 | small wording changes |
| html | 319c84c | tomkeaney | 2023-10-17 | Build site. |
| Rmd | 48c7e9f | tomkeaney | 2023-10-17 | Restarting the project |
| html | 92d2a26 | ausevo | 2023-03-08 | Build site. |
| Rmd | 5657ce7 | ausevo | 2023-03-08 | Better Table |
| html | 414b217 | ausevo | 2023-03-08 | Build site. |
| Rmd | 2501948 | ausevo | 2023-03-08 | Change to Tom’s analysis |
| html | c33ab4e | ausevo | 2023-03-08 | Build site. |
| Rmd | ca66d6a | ausevo | 2023-03-08 | Change to Tom’s analysis |
| html | daa935d | ausevo | 2023-03-08 | Build site. |
| Rmd | 78ece4f | ausevo | 2023-03-08 | Change to Tom’s analysis |
| html | 8f6bea5 | ausevo | 2023-03-03 | Build site. |
| Rmd | a382ebf | ausevo | 2023-03-03 | New version |
| html | b682a4f | ausevo | 2023-03-03 | Build site. |
| Rmd | 6577470 | ausevo | 2023-03-03 | New version |
| Rmd | 7bb248d | ausevo | 2023-01-31 | rework of line means results |
| html | 6d630f0 | ausevo | 2022-05-12 | Build site. |
| Rmd | 95f4a97 | ausevo | 2022-05-12 | including ideas section |
| html | 7bf0ee0 | ausevo | 2022-05-10 | Build site. |
| Rmd | 5b9883e | ausevo | 2022-05-10 | Add GCTA model outputs |
| Rmd | 8e2d8be | tkeaney | 2021-08-17 | Merge branch ‘master’ of https://github.com/tomkeaney/DGRP_sexual_conflict |
| Rmd | 4b095f0 | tkeaney | 2021-08-17 | edits |
| Rmd | 73f011c | lukeholman | 2021-08-06 | Tiny edit (digits=3) |
| html | 043902f | tkeaney | 2021-07-28 | Build site. |
| Rmd | 959123f | tkeaney | 2021-07-28 | trying to fix a bug |
| html | 769e48d | tkeaney | 2021-07-28 | Build site. |
| Rmd | 340c41b | tkeaney | 2021-07-28 | sexual dimorphism calculation added |
| html | 4018580 | tkeaney | 2021-07-28 | Build site. |
| Rmd | 65b2597 | tkeaney | 2021-07-28 | sexual dimorphism calculation added |
| html | 83a0c4f | tkeaney | 2021-07-28 | Build site. |
| Rmd | ddea55e | tkeaney | 2021-07-28 | sexual dimorphism calculation added |
| html | b364d30 | tkeaney | 2021-07-23 | Build site. |
| Rmd | e404cf9 | tkeaney | 2021-07-23 | progression on selection calculations |
| html | 32f107f | tkeaney | 2021-07-23 | Build site. |
| Rmd | 6aad8a8 | tkeaney | 2021-07-23 | progression on selection calculations |
| html | 11a8390 | tkeaney | 2021-07-07 | Build site. |
| Rmd | 42b8f12 | tkeaney | 2021-07-07 | Get the site up and running |
\(~\)
First load the packages and build helper functions
library(tidyverse) # for tidy coding
library(MetBrewer) # for many nice colour palettes
library(rcartocolor) # more cool colours
library(kableExtra) # for scrolling tables
library(DT) # for interactive tables
library(patchwork) # to join multiple plots nicely
library(brms) # for bayesian models
library(tidybayes) # for more bayesian things
library(bayestestR) # for the pd metric
library(broom) # convert results of functions into tables
library(ggtext) # for markdown features in ggplot
library(ggrepel) # for plot labels in ggplot
library(ggnewscale) # to reset scales in ggplot
library(bigsnpr)
# Create a function to build HTML searchable tables
my_data_table <- function(df){
datatable(
df, rownames=FALSE,
autoHideNavigation = TRUE,
extensions = c("Scroller", "Buttons"),
options = list(
autoWidth = TRUE,
dom = 'Bfrtip',
deferRender=TRUE,
scrollX=TRUE, scrollY=1000,
scrollCollapse=TRUE,
buttons =
list('pageLength', 'colvis', 'csv', list(
extend = 'pdf',
pageSize = 'A4',
orientation = 'landscape',
filename = 'Trait_data')),
pageLength = 100
)
)
}
\(~\)
We conducted a near-exhaustive search of the literature until January 2022, to obtain line mean estimates and associated meta-data for quantitative traits that have been measured in the DGRP. We did not include data collected for traits that had been measured in heterozygous combinations of multiple DGRP lines. In total, we identified 125 studies that reported line means or raw data for 2115 phenotypic traits. To ready the data for analysis, we grouped trait values by trait and sex and standardised the data to have mean = 0 and sd = 1. Sex-specific standardised line means for each trait were then combined with their standardised fitness estimates, obtained from Wong and Holman (2023). We also include some helpful metadata for downstream analysis. We then pruned the dataset to only include traits that have been measured in single-sex cohorts, in 80 or more lines. We also removed traits included in two large datasets on the microbiome and metabolome (Everett et al. 2020 and Jin et al. 2020).
# load in the data, note that traits have already been standardised
DGRP_data <-
left_join(
read_csv("data/input/all.dgrp.phenos_scaled.csv") %>%
mutate(line = as.factor(line)),
read_csv("data/input/meta_data_for_all_traits.csv") %>%
group_by(Reference) %>%
mutate(study_ID = as.factor(cur_group_id()),
Pooled = if_else(Sex == "Pooled", "Yes", "No"))) %>%
left_join(read_rds("data/input/trait_names.rds"))
# Apply the selection criteria with the filter() function, then add a column each for female and male fitness
trait_data <-
DGRP_data %>%
filter(!str_detect(Trait, "fitness"),
`# lines measured` >= 80 &
Pooled != "Yes" &
Reference != "Jin et al (2020) PLOS Genetics" &
Reference != "Everett et al (2020) Genome Research" &
!str_detect(Trait, "dopamine.response.to.paraquat.2021.m")) %>% # data for this trait was entered into the database incorrectly, it was only measured in 10 lines so should not be included in our analysis
# join the early life fitness data from Wong and Holman
left_join(
DGRP_data %>%
filter(str_detect(Trait, "fitness.early.life")) %>%
select(line, Trait, trait_value) %>%
pivot_wider(names_from = Trait, values_from = trait_value) %>%
rename(female_fitness = fitness.early.life.f, male_fitness = fitness.early.life.m))
clean_meta_data <-
trait_data %>%
select(Trait_nice, Trait, Life_stage, `Trait guild`, study_ID, Trait_nice, Reference, `Trait description`) %>%
distinct(Trait, .keep_all = TRUE) %>%
mutate(Trait_nice = case_when(Trait_nice == "DDT resistance mortality" ~ "DDT susceptibility",
.default = Trait_nice)) # make confusing name simpler
\(~\)
Table S1. Traits included for analysis
# how many studies did we start with?
table_data <- left_join(
DGRP_data %>%
distinct(Trait, .keep_all = TRUE) %>%
mutate(Measured_in = case_when(Sex == "Female" ~ "Females",
Sex == "Male" ~ "Males",
Sex == "Pooled" ~ "Mixed sex cohorts")) %>%
select(Trait_nice, Trait, Reference, `# lines measured`, Measured_in, `Trait description`),
trait_data %>%
distinct(Trait) %>%
mutate(Included = "Yes")
) %>%
mutate(Included = if_else(is.na(Included), "No", Included)) %>%
select(-Trait) %>%
rename(Trait = Trait_nice, `Measured in` = Measured_in) %>%
arrange(desc(Included))
# Create a function to build HTML searchable tables
my_data_table(table_data %>%
select(Trait, Reference, `# lines measured`, `Measured in`, Included, `Trait description`))
\(~\)
Find the number of traits and studies included in our analysis.
# how many studies and traits do we have after filtering?
num_unique_traits <- table_data %>% filter(Included == "Yes") %>% nrow()
# in females
num_unique_traits_f <- table_data %>% filter(Included == "Yes" & `Measured in` == "Females") %>% nrow()
# in males
num_unique_traits_m <- table_data %>% filter(Included == "Yes" & `Measured in` == "Males") %>% nrow()
# how many studies are they measured across in total?
num_unique_studies <- table_data %>% filter(Included == "Yes") %>% distinct(Reference) %>% nrow()
# in females
num_unique_studies_f <- table_data %>% filter(Included == "Yes" & `Measured in` == "Females") %>% distinct(Reference) %>% nrow()
# in males
num_unique_studies_m <- table_data %>% filter(Included == "Yes" & `Measured in` == "Males") %>% distinct(Reference) %>% nrow()
After this selection process, 474 remain, that were measured across 76 studies. There are 232 measured in females and 242 in males across 56 and 54 respectively.
\(~\)
\(~\)
We assume that line represents a genetically homogeneous, homozygous genotype, such that the line mean for a particular trait equals the breeding value for that trait for that genotype. We can therefore estimate the genetic variance for a given trait in the DGRP as the variance in line means (Mackay, 2012).
We can then estimate the response to selection (\(R\)) for a trait using Robertson’s Secondary Theorem of Natural Selection (also known as the Robertson covariance; Robertson, 1968), which states that \(R\) for a given trait is equivalent to the additive genetic covariance between the trait (\(A_z\)) and relative fitness (\(A_w\)), such that \(R = \sigma(A_w, Az)\). We can estimate this genetic covariance as the correlation in breeding values for traits with relative fitness across DGRP lines.
In species with two sexes and obligate sexual reproduction, half of the response to selection is due to selection on females and half to selection on males (e.g. Lande, 1980). We can therefore partition the secondary theorem into the response via each sex:
\[R = \frac{\sigma(A_w^F, A_z) + \sigma(A_w^M, A_z)}{2}\]
Furthermore, in quantitative genetics it is common to regard the same trait expressed in males or females (e.g. female body size and male body size) as separate traits. For convenience in the code below, we denote the response to selection of a female trait (\(R_F\)) or a male trait (\(R_M\)) as:
\[R_F = \frac{R_{FF} + R_{MF}}{2}\]
\[R_M = \frac{R_{MM} + R_{FM}}{2}\]
Where \(R_{FF}\) is the response of a female trait to selection on females, \(R_{MF}\) is the response of a female trait to selection on males, \(R_{MM}\) is the response of a male trait to selection on males, and \(R_{FM}\) is the response of a male trait to selection on females. These four quantities can be found via the partitioned Robertson covariance as follows:
\[R_{FF} = \sigma(A_w^F, A_z^F)\]
\[R_{MF} = \sigma(A_w^M, A_z^F)\]
\[R_{MM} = \sigma(A_w^M, A_z^M)\]
\[R_{FM} = \sigma(A_w^F, A_z^M)\]
Note that even though female traits are not expressed in males and therefore are not under selection in males, selection on males can still generate a response to selection in female traits if there is genetic covariance between these female traits and male traits that are under selection in males (and similarly, male traits can respond to selection on females due to between-sex genetic correlations). We therefore expect \(R_{MF}\) and \(R_{FM}\) to commonly be non-zero.
\(~\)
To estimate the covariance between \(A_w\) and \(A_z\), we fitted bivariate Bayesian linear
models using the brms package (Bürkner, 2017) for
R version 4.2.2. For each combination of trait and sex, we
used line means for the focal trait and the fitness of the focal sex as
the two response variables and fitted an intercept-only Gaussian model.
Each model returned a posterior distribution of the residual correlation
between trait and fitness, which for data expressed in standard units is
equivalent to the covariance.
Build functions to run the models
# RFF estimates
female_traits <- trait_data %>% filter(Sex == "Female")
trait_list_female <- unique(female_traits$Trait) # an input to the map_dfr() function that we'll need in a few chunks time
# code the model structure we will use for all traits using one example - `flight.performance.f`. We can then use the update() function to run this model many times, once for each trait measured in females. update() makes this process many times faster, because the model can immediately start sampling, without the need to recompile.
RFF_model <-
brm(data = female_traits %>% filter(Trait == "flight.performance.f"),
family = gaussian,
bf(mvbind(female_fitness, trait_value) ~ 1) + set_rescor(TRUE),
prior = c(prior(normal(0, 0.1), class = Intercept, resp = femalefitness),
prior(normal(0, 0.1), class = Intercept, resp = traitvalue),
prior(normal(1, 0.1), class = sigma, resp = femalefitness),
prior(normal(1, 0.1), class = sigma, resp = traitvalue),
prior(lkj(2), class = rescor)),
chains = 4, cores = 4, iter = 6000, warmup = 2000,
seed = 1, file = "fits/RFF_test_model")
# make a function to update the model and the posterior sample output with the 'selected trait'
RFF_calculator <- function(selected_trait){
data <- female_traits %>% filter(Trait == selected_trait)
model <- update(
RFF_model, newdata = data,
chains = 4, cores = 4, iter = 6000, warmup = 2000,
seed = 1)
posterior <-
as_draws_df(model) %>%
rename(Response_to_selection_female = rescor__femalefitness__traitvalue) %>%
mutate(Trait = selected_trait) %>%
select(Trait, Response_to_selection_female) %>%
as_tibble()
posterior
}
# RMF estimates
RMF_model <-
brm(data = female_traits %>% filter(Trait == "flight.performance.f"),
family = gaussian,
bf(mvbind(male_fitness, trait_value) ~ 1) + set_rescor(TRUE),
prior = c(prior(normal(0, 0.1), class = Intercept, resp = malefitness),
prior(normal(0, 0.1), class = Intercept, resp = traitvalue),
prior(normal(1, 0.1), class = sigma, resp = malefitness),
prior(normal(1, 0.1), class = sigma, resp = traitvalue),
prior(lkj(2), class = rescor)),
chains = 4, cores = 4, iter = 6000, warmup = 2000,
seed = 1, file = "fits/RMF_test_model")
# make a function to update the model and the posterior sample output with your desired trait
RMF_calculator <- function(selected_trait){
data <- female_traits %>% filter(Trait == selected_trait)
model <- update(
RMF_model, newdata = data,
chains = 4, cores = 4, iter = 6000, warmup = 2000,
seed = 1)
posterior <-
as_draws_df(model) %>%
rename(Response_to_selection_male = rescor__malefitness__traitvalue) %>%
mutate(Trait = selected_trait) %>%
select(Trait, Response_to_selection_male) %>%
as_tibble()
posterior
}
\(~\)
\(~\)
# RMM estimates
male_traits <- trait_data %>% filter(Sex == "Male")
trait_list_male <- unique(male_traits$Trait)
RMM_model <-
brm(data = male_traits %>% filter(Trait == "flight.performance.m"),
family = gaussian,
bf(mvbind(male_fitness, trait_value) ~ 1) + set_rescor(TRUE),
prior = c(prior(normal(0, 0.1), class = Intercept, resp = malefitness),
prior(normal(0, 0.1), class = Intercept, resp = traitvalue),
prior(normal(1, 0.1), class = sigma, resp = malefitness),
prior(normal(1, 0.1), class = sigma, resp = traitvalue),
prior(lkj(2), class = rescor)),
chains = 4, cores = 4, iter = 6000, warmup = 2000,
seed = 1, file = "fits/RMM_test_model")
# make a function to update the model and the posterior sample output with your desired trait
RMM_calculator <- function(selected_trait){
data <- male_traits %>% filter(Trait == selected_trait)
model <- update(
RMM_model, newdata = data,
chains = 4, cores = 4, iter = 6000, warmup = 2000,
seed = 1)
posterior <-
as_draws_df(model) %>%
rename(Response_to_selection_male = rescor__malefitness__traitvalue) %>%
mutate(Trait = selected_trait) %>%
select(Trait, Response_to_selection_male) %>%
as_tibble()
posterior
}
# RFM estimates
RFM_model <-
brm(data = male_traits %>% filter(Trait == "flight.performance.m"),
family = gaussian,
bf(mvbind(female_fitness, trait_value) ~ 1) + set_rescor(TRUE),
prior = c(prior(normal(0, 0.1), class = Intercept, resp = femalefitness),
prior(normal(0, 0.1), class = Intercept, resp = traitvalue),
prior(normal(1, 0.1), class = sigma, resp = femalefitness),
prior(normal(1, 0.1), class = sigma, resp = traitvalue),
prior(lkj(2), class = rescor)),
chains = 4, cores = 4, iter = 6000, warmup = 2000,
seed = 1, file = "fits/RFM_test_model")
# make a function to update the model and the posterior sample output with your desired trait
RFM_calculator <- function(selected_trait){
data <- male_traits %>% filter(Trait == selected_trait)
model <- update(
RFM_model, newdata = data,
chains = 4, cores = 4, iter = 6000, warmup = 2000,
seed = 1)
posterior <-
as_draws_df(model) %>%
rename(Response_to_selection_female = rescor__femalefitness__traitvalue) %>%
mutate(Trait = selected_trait) %>%
select(Trait, Response_to_selection_female) %>%
as_tibble()
posterior
}
\(~\)
Run the models using RFF_calculator,
RMF_calculator, RMM_calculator and
RFM_calculator
# run the RFF function
Run_function <- FALSE # change this to TRUE to run the models
if(Run_function){
RFF <- map_dfr(trait_list_female, RFF_calculator) # map_dfr returns a data frame created by row-binding each output
write_csv(RFF, file = "data/organismal_level_output/RFF.csv")
} else RFF <- read_csv("data/organismal_level_output/RFF.csv")
# run the RMF function
if(Run_function){
RMF <- map_dfr(trait_list_female, RMF_calculator)
write_csv(RMF, file = "data/organismal_level_output/RMF.csv")
} else RMF <- read_csv("data/organismal_level_output/RMF.csv")
# run the RMM function
if(Run_function){
RMM <- map_dfr(trait_list_male, RMM_calculator)
write_csv(RMM, file = "data/organismal_level_output/RMM.csv")
} else RMM <- read_csv("data/organismal_level_output/RMM.csv")
# run the RFM function
if(Run_function){
RFM <- map_dfr(trait_list_male, RFM_calculator)
write_csv(RFM, file = "data/organismal_level_output/RFM.csv")
} else RFM <- read_csv("data/organismal_level_output/RFM.csv")
\(~\)
R_female_traits <-
bind_cols(RFF, RMF) %>%
rename(Trait = Trait...1) %>%
mutate(Trait_Sex = "Female") %>%
select(-(Trait...3))
p1 <-
R_female_traits %>%
group_by(Trait) %>%
mutate(avg_R = median(Response_to_selection_female)) %>%
ggplot(aes(Response_to_selection_female, fct_reorder(Trait, avg_R))) +
stat_interval(.width = c(0.05, 0.66, 0.95),
height = 1, show.legend = F) +
rcartocolor::scale_color_carto_d(palette = "Purp") +
coord_cartesian(xlim = c(-0.5, 0.5)) +
geom_vline(linetype = 2, xintercept = 0, linewidth = 1) +
labs(x = "Response to selection in females",
y = "Trait expressed in females") +
theme_bw() +
theme(legend.position = "none",
panel.grid.minor = element_blank(),
text = element_text(size=14),
axis.text.y = element_text(size = 8))
p2 <-
R_female_traits %>%
group_by(Trait) %>%
mutate(avg_R = median(Response_to_selection_male)) %>%
ggplot(aes(Response_to_selection_male, fct_reorder(Trait, avg_R))) +
stat_interval(.width = c(0.05, 0.66, 0.95),
height = 1, show.legend = F) +
rcartocolor::scale_color_carto_d(palette = "Peach") +
coord_cartesian(xlim = c(-0.5, 0.5)) +
geom_vline(linetype = 2, xintercept = 0, linewidth = 1) +
labs(x = "Response to selection in males",
y = "Trait expressed in females") +
theme_bw() +
theme(legend.position = "none",
panel.grid.minor = element_blank(),
text = element_text(size=14),
axis.text.y = element_text(size = 8))
p1 + p2 +
plot_annotation(tag_levels = 'a')

Figure S1. The estimated response to selection in a females and b males for all traits measured in females. Innermost bands approximate the median, while outer bands show the 66 and 95% credible intervals.
\(~\)
R_male_traits <-
bind_cols(RFM, RMM) %>%
rename(Trait = Trait...1) %>%
mutate(Trait_Sex = "Male") %>%
select(-(Trait...3))
p3 <-
R_male_traits %>%
group_by(Trait) %>%
mutate(avg_R = median(Response_to_selection_male)) %>%
ggplot(aes(Response_to_selection_male, fct_reorder(Trait, avg_R))) +
stat_interval(.width = c(0.05, 0.66, 0.95),
height = 1, show.legend = F) +
rcartocolor::scale_color_carto_d(palette = "Peach") +
coord_cartesian(xlim = c(-0.5, 0.5)) +
geom_vline(linetype = 2, xintercept = 0, linewidth = 1) +
labs(x = "Response to selection in males",
y = "Trait expressed in males") +
theme_bw() +
theme(legend.position = "none",
panel.grid.minor = element_blank(),
text = element_text(size=14),
axis.text.y = element_text(size = 8))
p4 <-
R_male_traits %>%
group_by(Trait) %>%
mutate(avg_R = median(Response_to_selection_female)) %>%
ggplot(aes(Response_to_selection_female, fct_reorder(Trait, avg_R))) +
stat_interval(.width = c(0.05, 0.66, 0.95),
height = 1, show.legend = F) +
rcartocolor::scale_color_carto_d(palette = "Purp") +
coord_cartesian(xlim = c(-0.5, 0.5)) +
geom_vline(linetype = 2, xintercept = 0, linewidth = 1) +
labs(x = "Response to selection in females",
y = "Trait expressed in males") +
theme_bw() +
theme(legend.position = "none",
panel.grid.minor = element_blank(),
text = element_text(size=14),
axis.text.y = element_text(size = 8))
p3 + p4 +
plot_annotation(tag_levels = 'a')

Figure S2. The response to selection in a males and b females for all traits measured in males. Innermost bands approximate the median, while outer bands show the 66 and 95% credible intervals.
\(~\)
Combine the data frames and estimate the overall expected response to selection as
\[R = \frac{R_{F} + R_{M}}{2}\]
In the same code chunk, for each trait, we also calculate the difference in the response to selection acting on females and males as \(\Delta R = R_{FF} - R_{MF}\) for female traits or \(\Delta R = R_{FM} - R_{MM}\) for male traits. This difference facilitates identification of traits where the response to selection is very different between sexes (in magnitude, and possibly also in sign). Note that \(\Delta R\) alone does not reveal whether a trait is sexually antagonistic, because sexually concordant selection that varies in strength between sexes (or sex-limited selection) can result in a high value of \(\Delta R\).
R_female_traits <-
R_female_traits %>%
mutate(R_overall = (Response_to_selection_female + Response_to_selection_male)/2,
R_diff = Response_to_selection_female - Response_to_selection_male) %>%
select(Trait, Trait_Sex, everything())
R_male_traits <-
R_male_traits %>%
mutate(R_overall = (Response_to_selection_female + Response_to_selection_male)/2,
R_diff = Response_to_selection_female - Response_to_selection_male) %>%
select(Trait, Trait_Sex, everything()) %>%
filter(Trait != "dopamine.response.to.paraquat.2021.m")
R_all_traits <- bind_rows(R_female_traits, R_male_traits)
R_long_form <-
R_all_traits %>%
select(1:4) %>%
pivot_longer(cols = 3:4, names_to = "Fitness_Sex", values_to = "R_metric") %>%
mutate(Fitness_Sex = case_when(Fitness_Sex == "Response_to_selection_female" ~ "Female",
Fitness_Sex == "Response_to_selection_male" ~ "Male"))
\(~\)
We use the brms hypothesis() function to
compute evidence ratios, where our hypothesis is \(R\) (and its various related metrics) != 0.
When the hypothesis is one-sided, this is the posterior probability
under the hypothesis against its alternative. That is, when the
hypothesis is of the form R > 0, the evidence ratio is the ratio of
the posterior probability of R > 0 and the posterior probability of R
< 0. In this example, values greater than one indicate that the
evidence in favour of R > 0 is larger than evidence in favour of R
< 0. In contrast, values smaller than one indicate that there is
greater evidence in favour of R < 0 than R > 0. More on the
hypothesis function can be found here.
We consider evidence ratios > 39 or < 1/39 as biologically notable
(without accounting for multiple testing), as our hypothesis is
two-directional (negative and positive values of \(R\) are of interest to us) and these values
closely correspond to the familiar frequentist p-value = 0.05 (Makowski
et al, 2019).
# build a function to calculate the evidence ratio, get the relevant info from the output and convert it to a tibble
find_evidence_ratios <- function(Trait_name, specify_hypothesis){
x <- hypothesis(R_all_traits %>%
filter(Trait == Trait_name),
specify_hypothesis)
x <- x$hypothesis
x %>% as_tibble() %>%
select(Evid.Ratio) %>%
mutate(Trait = Trait_name)
}
# list of traits we need evidence ratios for.
all_traits_list <- R_all_traits %>% distinct(Trait) %>% pull(Trait)
# calculate evidence ratios
if(!file.exists("data/organismal_level_output/organismal_level_evidence_ratios.csv")){
R_evidence_ratios_female <-
map_dfr(all_traits_list,
find_evidence_ratios,
"Response_to_selection_female > 0") %>%
rename(Female_Evid_Ratio = Evid.Ratio)
R_evidence_ratios_male <-
map_dfr(all_traits_list,
find_evidence_ratios,
"Response_to_selection_male > 0") %>%
rename(Male_Evid_Ratio = Evid.Ratio)
R_evidence_ratios_overall <-
map_dfr(all_traits_list,
find_evidence_ratios,
"R_overall > 0") %>%
rename(Overall_Evid_Ratio = Evid.Ratio)
R_evidence_ratios_diff <-
map_dfr(all_traits_list,
find_evidence_ratios,
"R_diff > 0") %>%
rename(Diff_Evid_Ratio = Evid.Ratio)
}
Next, we manually calculate evidence ratios for sexually concordant responses to selection by:
Using the p_direction function from the
bayestestR package, we find the proportion of the posterior
distribution (the posterior probability) that is of the median’s sign
for each trait in each sex.
Using the p_direction output, we calculate the
probability that a trait has positive \(R\), \(P(pos)\) or negative \(R\), \(P(neg) = 1
- P(pos)\).
Find \(P(concord) = P(pos)_f P(pos)_m + P(neg)_f P(neg)_m\)
Find \(P(antag) = P(pos)_f P(neg)_m + P(neg)_f P(pos)_m\)
Take the ratio of \(P(concord)\) and \(P(antag)\)
Evidence ratios > 1 indicate that the response to selection is more likely to be sexually concordant, whereas ratios < 1 indicate a sexually antagonistic response has higher probability. As indicated by step 3, a concordant response is possible when the trait responds to selection in either a negative or positive direction in both sexes. Similarly, step 4 shows that an antagonistic response is also possible via two paths. Hence, we once again consider evidence ratios > 39 or < 1/39 as biologically notable (without accounting for multiple testing).
pd_function <-function(Trait_name){
R_all_traits %>% filter(Trait == Trait_name) %>%
select(Trait, Response_to_selection_female, Response_to_selection_male) %>%
p_direction(Response_to_selection_female, method = "direct", null = 0) %>%
as_tibble() %>%
pivot_wider(names_from = Parameter, values_from = pd) %>%
mutate(Trait = Trait_name) %>%
rename(pd_female = Response_to_selection_female,
pd_male = Response_to_selection_male)
}
if(!file.exists("data/organismal_level_output/organismal_level_evidence_ratios.csv")){
pd_data <- map_dfr(all_traits_list, pd_function)
Trait_medians <-
R_all_traits %>%
group_by(Trait) %>%
summarise(median_female = median(Response_to_selection_female),
median_male = median(Response_to_selection_male)) %>%
ungroup()
evidence_ratios_concord <-
left_join(pd_data, Trait_medians) %>%
mutate(prob_pos_female = if_else(median_female > 0, pd_female, 1 - pd_female),
prob_pos_male = if_else(median_male > 0, pd_male, 1 - pd_male)) %>%
select(Trait, prob_pos_female, prob_pos_male) %>%
# Calculate the probabilities that beta_i and beta_j have the same/opposite signs
mutate(p_sex_concord = prob_pos_female * prob_pos_male +
(1 - prob_pos_female) * (1 - prob_pos_male),
p_sex_antag = prob_pos_female * (1 - prob_pos_male) +
(1 - prob_pos_female) * prob_pos_male) %>%
# Find the ratios of these two probabilities (i.e. the "evidence ratios")
mutate(Concord_Evid_Ratio = p_sex_concord / p_sex_antag)
}
Join the evidence ratios into a single tibble and save as a
.csv file for fast loading. Create new columns that
indicate whether the evidence ratio indicate biologically
notability.
if(!file.exists("data/organismal_level_output/organismal_level_evidence_ratios.csv")){
organismal_level_evidence_ratios <-
R_evidence_ratios_female %>%
left_join(R_evidence_ratios_male) %>%
left_join(R_evidence_ratios_overall) %>%
left_join(R_evidence_ratios_diff) %>%
left_join(evidence_ratios_concord %>% select(Trait, Concord_Evid_Ratio)) %>%
mutate(Female_notable = if_else(Female_Evid_Ratio > 39 | Female_Evid_Ratio < 1/39, "Yes", "No"),
Male_notable = if_else(Male_Evid_Ratio > 39 | Male_Evid_Ratio < 1/39, "Yes", "No"),
Overall_notable = if_else(Overall_Evid_Ratio > 39 | Overall_Evid_Ratio < 1/39, "Yes", "No"),
Diff_notable = if_else(Diff_Evid_Ratio > 39 | Diff_Evid_Ratio < 1/39, "Yes", "No"),
Concord_notable = if_else(Concord_Evid_Ratio > 39 | Concord_Evid_Ratio < 1/39, "Yes", "No")) %>%
left_join(clean_meta_data) %>%
left_join(R_all_traits %>% distinct(Trait, Trait_Sex)) %>%
mutate(Trait = fct_reorder(Trait, `Trait guild`)) %>%
arrange(`Trait guild`, Trait) %>%
mutate(position = 1:n(),
Trait_Sex = if_else(Trait_Sex == "Female", "Traits measured in females",
"Traits measured in males")) %>%
select(Trait, everything())
write_csv(evidence_ratios, "data/organismal_level_output/organismal_level_evidence_ratios.csv")
} else organismal_level_evidence_ratios <- read_csv("data/organismal_level_output/organismal_level_evidence_ratios.csv")
\(~\)
Build the evidence ratio plots
# work out where to place the x axis labels - we want them in the middle of their trait guild
x_labels <- organismal_level_evidence_ratios %>%
group_by(`Trait guild`) %>%
summarise(position = min(position) + (max(position) - min(position))/2)
# get colours for each trait guild
guild_pal <- c("#6e7cb9", "#7bbcd5", "#d0e2af", "#f5db99",
"#e89c81", "#d2848d", "#6e7cb9", "#7bbcd5",
"#d0e2af","#f5db99","#e89c81", "#d2848d",
"#6e7cb9", "#7bbcd5")
plot_a <-
organismal_level_evidence_ratios %>%
ggplot(aes(x = position, y = log2(Overall_Evid_Ratio))) +
geom_point(aes(colour = `Trait guild`),
size = 4.5, alpha = 1, shape = 17) +
geom_hline(yintercept = log2(1), linetype = 2, linewidth = 1) +
geom_hline(yintercept = log2(39), linetype = 3, linewidth = 0.9) +
geom_hline(yintercept = log2(1/39), linetype = 3, linewidth = 0.9) +
#scale_shape_manual(values = c(25, 24)) +
scale_colour_manual(values = guild_pal) +
scale_x_continuous(breaks = x_labels$position, labels = x_labels$`Trait guild`) +
scale_y_continuous(expand = c(0.01, 0), breaks = c(-12, -10, -8, -6, -4, -2, 0, 2, 4, 6, 8, 10, 12),
labels = c(paste("1/",2 ^ c(12, 10, 8, 6, 4, 2), sep = ""), 2 ^ c(0, 2,4,6,8,10,12)), limits = c(-14, 14)) +
labs(x = NULL,
#y = "R~F~ log2(ER)",
y = "_R_ log~2~(ER)") +
geom_label_repel(data = organismal_level_evidence_ratios %>%
filter(Overall_Evid_Ratio > 200 | Overall_Evid_Ratio < 1/200),
aes(label = Trait_nice),
fill = "white",
size = 3, alpha = 0.9,
min.segment.length = 0, seed = 1,
box.padding = 0.5, point.padding = 0.5,
max.overlaps = 30) +
theme_minimal() +
theme(legend.position = "none",
panel.grid.major.x = element_blank(),
panel.grid.minor.x = element_blank(),
panel.grid.minor.y = element_blank(),
axis.text.x=element_blank(),
axis.ticks.x=element_blank(),
axis.title.x = element_text(size = 15),
axis.title.y = element_markdown(size = 15))
plot_b <-
organismal_level_evidence_ratios %>%
ggplot(aes(x = position, y = log2(Female_Evid_Ratio))) +
geom_point(aes(colour = `Trait guild`),
size = 4.5, alpha = 1, shape = 17) +
geom_hline(yintercept = log2(1), linetype = 2, linewidth = 1) +
geom_hline(yintercept = log2(39), linetype = 3, linewidth = 0.9) +
geom_hline(yintercept = log2(1/39), linetype = 3, linewidth = 0.9) +
#scale_shape_manual(values = c(25, 24)) +
scale_colour_manual(values = guild_pal) +
scale_x_continuous(breaks = x_labels$position, labels = x_labels$`Trait guild`) +
scale_y_continuous(expand = c(0.01, 0), breaks = c(-12, -10, -8, -6, -4, -2, 0, 2, 4, 6, 8, 10, 12),
labels = c(paste("1/",2 ^ c(12, 10, 8, 6, 4, 2), sep = ""), 2 ^ c(0, 2,4,6,8,10,12)), limits = c(-14, 14)) +
labs(x = NULL,
#y = "R~F~ log2(ER)",
y = "_R_ on females log~2~(ER)") +
geom_label_repel(data = organismal_level_evidence_ratios %>%
filter(Female_Evid_Ratio > 100 | Female_Evid_Ratio < 1/100),
aes(label = Trait_nice),
fill = "white",
size = 3, alpha = 0.9,
min.segment.length = 0, seed = 1,
box.padding = 0.5, point.padding = 0.5,
max.overlaps = 30) +
theme_minimal() +
theme(legend.position = "none",
panel.grid.major.x = element_blank(),
panel.grid.minor.x = element_blank(),
panel.grid.minor.y = element_blank(),
axis.text.x=element_blank(),
axis.ticks.x=element_blank(),
axis.title.x = element_text(size = 15),
axis.title.y = element_markdown(size = 15))
plot_c <-
organismal_level_evidence_ratios %>%
ggplot(aes(x = position, y = log2(Male_Evid_Ratio))) +
geom_point(aes(colour = `Trait guild`),
size = 4.5, alpha = 1, shape = 17) +
geom_hline(yintercept = log2(1), linetype = 2, linewidth = 1) +
geom_hline(yintercept = log2(39), linetype = 3, linewidth = 0.9) +
geom_hline(yintercept = log2(1/39), linetype = 3, linewidth = 0.9) +
#scale_shape_manual(values = c(25, 24)) +
scale_colour_manual(values = guild_pal) +
scale_x_continuous(breaks = x_labels$position, labels = x_labels$`Trait guild`) +
scale_y_continuous(expand = c(0.01, 0), breaks = c(-12, -10, -8, -6, -4, -2, 0, 2, 4, 6, 8, 10, 12),
labels = c(paste("1/",2 ^ c(12, 10, 8, 6, 4, 2), sep = ""), 2 ^ c(0, 2,4,6,8,10,12)), limits = c(-14, 14)) +
labs(x = "Traits",
y = "_R_ on males log~2~(ER)") +
geom_label_repel(data = organismal_level_evidence_ratios %>%
filter(Male_Evid_Ratio > 100 | Male_Evid_Ratio < 1/100),
aes(label = Trait_nice),
fill = "white",
size = 3, alpha = 0.9,
min.segment.length = 0, seed = 1,
box.padding = 0.5, point.padding = 0.5,
max.overlaps = 30) +
theme_minimal() +
theme(legend.position = "none",
panel.grid.major.x = element_blank(),
panel.grid.minor.x = element_blank(),
panel.grid.minor.y = element_blank(),
axis.text.x=element_text(angle = 90,
vjust = 0.5,
size = 10),
axis.title.x = element_text(size = 15),
axis.title.y = element_markdown(size = 15))
(Figure_S3 <- plot_a / plot_b / plot_c +
plot_annotation(tag_levels = 'a'))

Figure S3. Triangles show evidence ratios (ERs) for a overall \(R\), b \(R\) on females (\(R_{FF}\) and \(R_{MF}\)) and c \(R\) on males (\(R_{FM}\) and \(R_{MM}\)). Traits are loosely organised into guilds to aid visualisation. The dashed lines indicate an evidence ratio of 1, where the probability that \(R\) > 0 is equal to the probability that \(R\) < 0. ER > 1 indicates a positive response is more likely, whereas ER < 1 indicates a negative response is more likely. The upper dotted line on each plot indicates an evidence ratio of 39, while the lower dotted line indicates an evidence ratio of 1/39; these correspond strongly with the frequentist \(p\) = 0.05. Traits with the most extreme evidence ratios are labelled.
plot_d <-
organismal_level_evidence_ratios %>%
ggplot(aes(x = position, y = log2(Diff_Evid_Ratio))) +
geom_point(aes(colour = `Trait guild`),
size = 4.5, alpha = 1, shape = 17) +
geom_hline(yintercept = log2(1), linetype = 2, linewidth = 1) +
geom_hline(yintercept = log2(39), linetype = 3, linewidth = 0.9) +
geom_hline(yintercept = log2(1/39), linetype = 3, linewidth = 0.9) +
scale_colour_manual(values = guild_pal) +
scale_x_continuous(breaks = x_labels$position, labels = x_labels$`Trait guild`) +
scale_y_continuous(expand = c(0.01, 0), breaks = c(-10, -8, -6, -4, -2, 0, 2, 4, 6, 8, 10),
labels = c(paste("1/",2 ^ c(10, 8, 6, 4, 2), sep = ""), 2 ^ c(0, 2,4,6,8,10)), limits = c(-12, 12)) +
labs(x = NULL,
y = "Δ_R_ log~2~(ER)") +
geom_label_repel(data = organismal_level_evidence_ratios %>%
mutate(Trait_nice = case_when(Trait_nice == "DDT resistance mortality" ~ "DDT susceptibility", .default = Trait_nice)) %>%
filter(Diff_Evid_Ratio > 39 | Diff_Evid_Ratio < 1/39),
aes(label = Trait_nice),
fill = "white",
size = 3, alpha = 0.9,
min.segment.length = 0, seed = 1,
box.padding = 0.5, point.padding = 0.5,
max.overlaps = 30) +
theme_minimal() +
theme(legend.position = "none",
panel.grid.major.x = element_blank(),
panel.grid.minor.x = element_blank(),
panel.grid.minor.y = element_blank(),
axis.text.x=element_blank(),
axis.ticks.x=element_blank(),
axis.title.x = element_text(size = 15),
axis.title.y = element_markdown(size = 15))
plot_e <-
organismal_level_evidence_ratios %>%
ggplot(aes(x = position, y = log2(Concord_Evid_Ratio))) +
geom_point(aes(colour = `Trait guild`),
size = 4.5, alpha = 1, shape = 17) +
geom_hline(yintercept = log2(1), linetype = 2, linewidth = 1) +
geom_hline(yintercept = log2(39), linetype = 3, linewidth = 0.9) +
geom_hline(yintercept = log2(1/39), linetype = 3, linewidth = 0.9) +
scale_colour_manual(values = guild_pal) +
scale_x_continuous(breaks = x_labels$position, labels = x_labels$`Trait guild`) +
scale_y_continuous(expand = c(0.01, 0), breaks = c(-10, -8, -6, -4, -2, 0, 2, 4, 6, 8, 10),
labels = c(paste("1/",2 ^ c(10, 8, 6, 4, 2), sep = ""), 2 ^ c(0, 2,4,6,8,10)), limits = c(-12, 12)) +
labs(x = "Trait",
y = "Sexual concordance log~2~(ER)") +
geom_label_repel(data = organismal_level_evidence_ratios %>%
filter(Concord_Evid_Ratio > 39 | Concord_Evid_Ratio < 1/39),
aes(label = Trait_nice),
fill = "white",
size = 3, alpha = 0.9,
min.segment.length = 0, seed = 1,
box.padding = 0.5, point.padding = 0.5,
max.overlaps = 30) +
theme_minimal() +
theme(legend.position = "none",
panel.grid.major.x = element_blank(),
panel.grid.minor.x = element_blank(),
panel.grid.minor.y = element_blank(),
axis.text.x=element_text(angle = 90,
vjust = 0.5,
size = 10),
axis.title.x = element_text(size = 15),
axis.title.y = element_markdown(size = 15))
(Figure_3 <- plot_d/ plot_e + plot_annotation(tag_levels = "a"))

Figure 3. Triangles show evidence ratios (ERs) for a \(\Delta R\) (\(R\) on females - \(R\) on males) and b sexual concordance, for all measured traits in the organismal level phenotype dataset. Traits are loosely organised into guilds to aid visualisation. The dashed lines indicate an evidence ratio of 1, where the probability that \(\Delta R\) > 0 is equal to the probability that \(\Delta R\) < 0, or that a trait is just as likely to have a sexually concordant response to selection (ER > 1) as a sexually antagonistic response (ER < 1). The upper dotted line on each plot indicates an evidence ratio of 39, while the lower dotted line indicates an evidence ratio of 1/39; these correspond strongly with the frequentist \(p\) = 0.05. Traits with the most extreme evidence ratios are labelled.
Table S2 / Supplementary dataset
S1: \(R\) partitions and
associated evidence ratios for the organismal level trait dataset. Trait
descriptions are provided in Table S1. Data can be downloaded as a
.csv file from the github
repository associated with this project.
all_ol_traits_summary <-
R_all_traits %>%
group_by(Trait, Trait_Sex) %>%
median_qi(Response_to_selection_female,
Response_to_selection_male,
R_overall,
R_diff) %>%
rename(`Sex trait was measured in` = Trait_Sex,
R = R_overall,
`R 2.5% CI` = R_overall.lower,
`R 97.5% CI` = R_overall.upper,
`R female` = Response_to_selection_female,
`R female 2.5% CI` = Response_to_selection_female.lower,
`R female 97.5% CI` = Response_to_selection_female.upper,
`R male` = Response_to_selection_male,
`R male 2.5% CI` = Response_to_selection_male.lower,
`R male 97.5% CI` = Response_to_selection_male.upper,
`Delta R` = R_diff,
`Delta R 2.5% CI` = R_diff.lower,
`Delta R 97.5% CI` = R_diff.upper) %>%
left_join(organismal_level_evidence_ratios %>% select(-contains("notable")) %>%
rename(`R ER` = Overall_Evid_Ratio, `R female ER` = Female_Evid_Ratio,
`R male ER` = Male_Evid_Ratio, `Delta R ER` = Diff_Evid_Ratio,
`Sexual concordance ER` = Concord_Evid_Ratio)) %>%
select(Trait, Trait_nice, `Sex trait was measured in`, R, `R 2.5% CI`, `R 97.5% CI`, `R ER`,
`R female`, `R female 2.5% CI`, `R female 97.5% CI`, `R female ER`,
`R male`, `R male 2.5% CI`, `R male 97.5% CI`, `R male ER`,
`Delta R`, `Delta R 2.5% CI`, `Delta R 97.5% CI`, `Delta R ER`, `Sexual concordance ER`) %>%
ungroup() %>%
mutate(across(c(4,5,6,8,9,10,12,13,14,16,17,18), ~round(.x, 3)),
across(ends_with("ER"), ~ round(.x, 4)))
write_csv(all_ol_traits_summary, "data/organismal_level_output/Supplementary_dataset_1.csv")
my_data_table(all_ol_traits_summary)
\(~\)
Fit the model to test whether \(|R|\) is affected by the sex fitness and trait values were measured in. Here, we take a meta-analysis type approach and weight each \(R\) estimate by the inverse of its standard error.
R_medians <-
R_long_form %>%
group_by(Trait, Fitness_Sex, Trait_Sex) %>%
summarise_draws("median", "sd", ~quantile(.x, probs = c(0.025, 0.975), na.rm = TRUE)) %>%
ungroup() %>%
select(-variable) %>%
left_join(clean_meta_data) %>%
mutate(absolute_R = abs(median),
Fitness_Sex = as.factor(Fitness_Sex),
Trait_Sex = as.factor(Trait_Sex),
Trait = as.factor(Trait))
median_R_model <-
brm(absolute_R | weights(1/sd) ~ 1 + Fitness_Sex * Trait_Sex + (1|study_ID) + (1|Trait),
family = brmsfamily(family = "Gamma"), # gamma is appropriate for the half-normal distribution created by taking the absolute
data = R_medians,
prior = c(prior(normal(-2.2, 1), class = Intercept),
prior(exponential(1), class = sd),
prior(exponential(1), class = shape),
prior(normal(0, 0.25), class = b)),
warmup = 4000, iter = 12000,
seed = 1, cores = 4, chains = 4,
control = list(adapt_delta = 0.9, max_treedepth = 10),
file = "fits/median_R_model")
print(median_R_model)
Family: gamma
Links: mu = log; shape = identity
Formula: absolute_R | weights(1/sd) ~ 1 + Fitness_Sex * Trait_Sex + (1 | study_ID) + (1 | Trait)
Data: R_medians (Number of observations: 950)
Draws: 4 chains, each with iter = 12000; warmup = 4000; thin = 1;
total post-warmup draws = 32000
Group-Level Effects:
~study_ID (Number of levels: 76)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept) 0.32 0.05 0.23 0.43 1.00 6142 10131
~Trait (Number of levels: 475)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept) 0.59 0.02 0.54 0.63 1.00 8421 12837
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Rhat
Intercept -2.70 0.06 -2.82 -2.58 1.00
Fitness_SexMale 0.18 0.02 0.14 0.22 1.00
Trait_SexMale -0.15 0.06 -0.27 -0.03 1.00
Fitness_SexMale:Trait_SexMale 0.11 0.03 0.05 0.17 1.00
Bulk_ESS Tail_ESS
Intercept 16404 18313
Fitness_SexMale 36862 23998
Trait_SexMale 11981 15341
Fitness_SexMale:Trait_SexMale 34545 23661
Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
shape 2.42 0.03 2.36 2.49 1.00 39872 21791
Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).
To check how does the gamma distribution fares at capturing the shape of the data, we use the model to recapitulate the data it was trained on. Given the shape of the predicted data is similar to the real data, there’s a high chance we’ve used an appropriate distribution family.
pp_check(median_R_model)

| Version | Author | Date |
|---|---|---|
| ad196cf | tomkeaney | 2024-04-11 |
\(~\)
Get model predictions and plot
new_data <- expand_grid(Fitness_Sex = c("Female", "Male"),
Trait_Sex = c("Female", "Male"))
R_fitted <-
fitted(median_R_model, newdata = new_data, re_formula = NA, summary = F) %>%
as.data.frame() %>%
rename(FemaleFitness_FemaleTrait = V1, FemaleFitness_MaleTrait = V2,
MaleFitness_FemaleTrait = V3, MaleFitness_MaleTrait = V4) %>%
as_tibble() %>%
mutate(percent_diff_female = ((MaleFitness_FemaleTrait - FemaleFitness_FemaleTrait) / FemaleFitness_FemaleTrait)*100,
percent_diff_male = ((MaleFitness_MaleTrait - FemaleFitness_MaleTrait)/ FemaleFitness_MaleTrait)*100) %>%
pivot_longer(cols = everything(), names_to = "Parameter", values_to = "R_mean") %>%
mutate(Fitness_Sex = case_when(str_detect(Parameter, "FemaleFitness") ~ "Female",
str_detect(Parameter, "MaleFitness") ~ "Male"),
Trait_Sex = case_when(str_detect(Parameter, "FemaleTrait") ~ "Female",
str_detect(Parameter, "MaleTrait") ~ "Male",
str_detect(Parameter, "percent_diff_female") ~ "Female",
str_detect(Parameter, "percent_diff_male") ~ "Male"))
R_p1 <-
R_fitted %>%
filter(!str_detect(Parameter, "percent")) %>%
ggplot(aes(x = R_mean, y = Trait_Sex, fill = Fitness_Sex)) + #fill = Fitness_Sex)) +
stat_slab(alpha = 0.9, shape = 21) +#,
labs(x = expression("Mean |"* italic(R) * "| for organismal level traits"),
y = "Phenotyped sex") +
scale_fill_manual(values = c(carto_pal(7, "Purp")[5], carto_pal(7, "Peach")[5])) +
coord_cartesian(xlim = c(0.045, 0.10)) +
theme_minimal() +
theme(panel.background = element_rect(fill='transparent'),
#panel.grid.major.y = element_blank(),
panel.grid.minor.y = element_blank(),
plot.background = element_rect(fill='transparent', color=NA),
legend.position = "none",
text = element_text(size=13))
R_p2 <-
R_fitted %>%
filter(str_detect(Parameter, "percent")) %>%
ggplot(aes(x = R_mean, y = Trait_Sex, fill = after_stat(x > 0))) +
stat_halfeye(.width = c(0.66, 0.95), alpha = 0.9, fill = carto_pal(7, "Emrld")[2],
point_interval = "median_qi", point_fill = "white",
shape = 21, point_size = 4, stroke = 1.2) +
geom_vline(xintercept = 0, linetype = 2, linewidth = 1.2) +
#coord_cartesian(xlim = c(-5, 50)) +
xlab(expression("% sex difference in |"*italic(R)* "|")) +
ylab("Phenotyped sex") +
coord_cartesian(xlim = c(-5, 45)) +
theme_minimal() +
theme(panel.background = element_rect(fill='transparent'),
#panel.grid.major.y = element_blank(),
panel.grid.minor.y = element_blank(),
plot.background = element_rect(fill='transparent', color=NA),
legend.position = "none",
text = element_text(size=13))
R_p1 + R_p2 +
plot_annotation(tag_levels = 'a')

\(~\)
Calculating stats for the results section
R_fitted %>%
group_by(Parameter) %>%
median_qi(R_mean)
# A tibble: 6 × 7
Parameter R_mean .lower .upper .width .point .interval
<chr> <dbl> <dbl> <dbl> <dbl> <chr> <chr>
1 FemaleFitness_FemaleTrait 0.0673 0.0598 0.0761 0.95 median qi
2 FemaleFitness_MaleTrait 0.0579 0.0513 0.0653 0.95 median qi
3 MaleFitness_FemaleTrait 0.0808 0.0716 0.0915 0.95 median qi
4 MaleFitness_MaleTrait 0.0776 0.0689 0.0876 0.95 median qi
5 percent_diff_female 20.1 15.3 25.1 0.95 median qi
6 percent_diff_male 34.1 28.7 39.7 0.95 median qi
\(~\)
We use the data archived on the dgrp2 website, which was provided by Huang et al. (2015). These data are levels of expression of genes measured across 185 DGRP lines. Huang et al’s data contains Y-linked genes that have higher/equal expression in females in all lines, presumably microarray issues/errors. To be conservative, we restrict our analyses to genes that are known to be on chromosomes that are present in both sexes. After data cleaning, we retain 14,268 genes in our analysis.
We also load gene annotation data using the org.Dm.eg.db
R package provided by BiocManager. The code used to produce
annotations is provided in the code subdirectory, in the
get_gene_annotations.R file.
Huang et al. 2015 replicated each gene expression measurement twice. To find line means, we simply take the average value for each line. We then standardise the expression of each gene to have \(\mu = 0\) and \(\sigma = 1\) (as with the ‘phenotypic’ traits.
# load in the gene names, we'll use these for plots and tables
Dmel_gene_names <- read_csv("data/input/all_dmel_genes.csv")
gene_annotations <- read_csv("data/input/gene_anntotations.csv")
# Helper to load all the Huang et al. expression data into tidy format
load_expression_data <- function(gene_annotations){
# Note: Huang et al's data contains Y-linked genes that have
# higher/equal expression in *females* in all lines, presumably microarray issues/errors.
# To be conservative, we restrict our analyses to genes that are known to be on a
# chromosomes that is present in both sexes
genes_allowed <- gene_annotations %>%
filter(chromosome %in% c("2L", "2R", "3L", "3R", "4", "X")) %>%
pull(FBID)
females <- read_delim("data/input/huang_transcriptome/dgrp.array.exp.female.txt", delim = " ") %>%
filter(gene %in% genes_allowed) %>%
gather(sample, expression, -gene) %>%
mutate(line = map_chr(str_extract_all(sample, "line_[:digit:]*"), ~ .x[1]),
replicate = map_chr(str_split(sample, ":"), ~ .x[2]),
sex = "Female")
males <- read_delim("data/input/huang_transcriptome/dgrp.array.exp.male.txt", delim = " ") %>%
filter(gene %in% genes_allowed) %>%
gather(sample, expression, -gene) %>%
mutate(line = map_chr(str_extract_all(sample, "line_[:digit:]*"), ~ .x[1]),
replicate = map_chr(str_split(sample, ":"), ~ .x[2]),
sex = "Male")
bind_rows(females, males) %>%
select(sample, line, sex, replicate, gene, expression)
}
expression_line_means <-
load_expression_data(gene_annotations) %>% # use the custom function to load expression data
mutate(line = str_remove(line, "line_"),
line = as.factor(line)) %>%
group_by(line, sex, gene) %>%
mutate(trait_value = mean(expression)) %>% # compute the average between replicates for each gene
ungroup() %>%
distinct(line, sex, gene, trait_value) %>%
group_by(sex, gene) %>% # scale the traits, specific to gene and sex
mutate(trait_value = as.numeric(scale(trait_value))) %>%
ungroup() %>%
rename(Sex = sex) %>%
# join the fitness data
left_join(trait_data %>% distinct(line, female_fitness, male_fitness))
# the transcriptome is large; memory is therefore a consistent problem with this analysis - it helps to clear often
#
invisible(gc())
\(~\)
\(~\)
To estimate the covariance between \(A_w\) and \(A_z\) (which in this case is a transcript),
we fitted bivariate Bayesian linear models using the brms
package (Bürkner, 2017) for R version 4.2.2. For each
combination of trait/transcript and sex, we used line means for the
focal trait/transcript and the fitness of the focal sex as the two
response variables and fitted an intercept-only Gaussian model. Each
model returned a posterior distribution of the residual correlation
between trait/transcript and fitness, which for data expressed in
standard units is equivalent to the covariance.
Build functions to run the models
# RFF estimates
female_transcripts <- expression_line_means %>% filter(Sex == "Female")
transcript_list <- unique(female_transcripts$gene) # an input to the map_dfr() function that we'll need in a few chunks time
# code the model structure we will use for all traits/transcripts using one example - `FBgn0000014`. We can then use the update() function to run this model many times, once for each trait/transcript measured in females. update() makes this process many times faster, because the model can immediately start sampling, without the need to recompile.
RFF_transcriptome_model <-
brm(data = female_transcripts %>% filter(gene == "FBgn0000014"),
family = gaussian,
bf(mvbind(female_fitness, trait_value) ~ 1) + set_rescor(TRUE),
prior = c(prior(normal(0, 0.1), class = Intercept, resp = femalefitness),
prior(normal(0, 0.1), class = Intercept, resp = traitvalue),
prior(normal(1, 0.1), class = sigma, resp = femalefitness),
prior(normal(1, 0.1), class = sigma, resp = traitvalue),
prior(lkj(2), class = rescor)),
chains = 4, cores = 4, iter = 4000, warmup = 2000,
seed = 1, file = "fits/RFF_transcriptome_model",
backend = "cmdstanr", refresh = 400)
# make a function to update the model and the posterior sample output with the 'selected trait'
RFF_transcriptome_calculator <- function(selected_gene){
data <- female_transcripts %>% filter(gene == selected_gene)
model <- update(
RFF_transcriptome_model, newdata = data,
chains = 4, cores = 4, iter = 4000, warmup = 2000,
seed = 1, backend = "cmdstanr", refresh = 400)
posterior <-
as_draws_df(model) %>%
rename(Response_to_selection_female = rescor__femalefitness__traitvalue) %>%
mutate(Trait = selected_gene) %>%
select(Trait, Response_to_selection_female) %>%
as_tibble()
posterior
}
# RMF estimates
RMF_transcriptome_model <-
brm(data = female_transcripts %>% filter(gene == "FBgn0000014"),
family = gaussian,
bf(mvbind(male_fitness, trait_value) ~ 1) + set_rescor(TRUE),
prior = c(prior(normal(0, 0.1), class = Intercept, resp = malefitness),
prior(normal(0, 0.1), class = Intercept, resp = traitvalue),
prior(normal(1, 0.1), class = sigma, resp = malefitness),
prior(normal(1, 0.1), class = sigma, resp = traitvalue),
prior(lkj(2), class = rescor)),
chains = 4, cores = 4, iter = 4000, warmup = 2000,
seed = 1, file = "fits/RMF_transcriptome_model",
backend = "cmdstanr", refresh = 400)
# make a function to update the model and the posterior sample output with your desired trait
RMF_transcriptome_calculator <- function(selected_gene){
data <- female_transcripts %>% filter(gene == selected_gene)
model <- update(
RMF_transcriptome_model, newdata = data,
chains = 4, cores = 4, iter = 4000, warmup = 2000,
seed = 1, backend = "cmdstanr", refresh = 400)
posterior <-
as_draws_df(model) %>%
rename(Response_to_selection_male = rescor__malefitness__traitvalue) %>%
mutate(Trait = selected_gene) %>%
select(Trait, Response_to_selection_male) %>%
as_tibble()
posterior
}
\(~\)
\(~\)
# RMM estimates
male_transcripts <- expression_line_means %>% filter(Sex == "Male")
RMM_transcriptome_model <-
brm(data = male_transcripts %>% filter(gene == "FBgn0000014"),
family = gaussian,
bf(mvbind(male_fitness, trait_value) ~ 1) + set_rescor(TRUE),
prior = c(prior(normal(0, 0.1), class = Intercept, resp = malefitness),
prior(normal(0, 0.1), class = Intercept, resp = traitvalue),
prior(normal(1, 0.1), class = sigma, resp = malefitness),
prior(normal(1, 0.1), class = sigma, resp = traitvalue),
prior(lkj(2), class = rescor)),
chains = 4, cores = 4, iter = 4000, warmup = 2000,
seed = 1, file = "fits/RMM_transcriptome_model",
backend = "cmdstanr", refresh = 400)
# make a function to update the model and the posterior sample output with your desired trait
RMM_transcriptome_calculator <- function(selected_gene){
data <- male_transcripts %>% filter(gene == selected_gene)
model <- update(
RMM_transcriptome_model, newdata = data,
chains = 4, cores = 4, iter = 4000, warmup = 2000,
seed = 1)
posterior <-
as_draws_df(model) %>%
rename(Response_to_selection_male = rescor__malefitness__traitvalue) %>%
mutate(Trait = selected_gene) %>%
select(Trait, Response_to_selection_male) %>%
as_tibble()
posterior
}
# RFM estimates
RFM_transcriptome_model <-
brm(data = male_transcripts %>% filter(gene == "FBgn0000014"),
family = gaussian,
bf(mvbind(female_fitness, trait_value) ~ 1) + set_rescor(TRUE),
prior = c(prior(normal(0, 0.1), class = Intercept, resp = femalefitness),
prior(normal(0, 0.1), class = Intercept, resp = traitvalue),
prior(normal(1, 0.1), class = sigma, resp = femalefitness),
prior(normal(1, 0.1), class = sigma, resp = traitvalue),
prior(lkj(2), class = rescor)),
chains = 4, cores = 4, iter = 4000, warmup = 2000,
seed = 1, file = "fits/RFM_transcriptome_model",
backend = "cmdstanr", refresh = 400)
# make a function to update the model and the posterior sample output with your desired trait
RFM_transcriptome_calculator <- function(selected_trait){
data <- male_transcripts %>% filter(gene == selected_trait)
model <- update(
RFM_transcriptome_model, newdata = data,
chains = 4, cores = 4, iter = 4000, warmup = 2000,
seed = 1)
posterior <-
as_draws_df(model) %>%
rename(Response_to_selection_female = rescor__femalefitness__traitvalue) %>%
mutate(Trait = selected_gene) %>%
select(Trait, Response_to_selection_female) %>%
as_tibble()
posterior
}
\(~\)
Run the models using RFF_transcriptome_calculator,
RMF_transcriptome_calculator,
RMM_transcriptome_calculator and
RFM_transcriptome_calculator.
This takes a fair bit of memory, so it might be necessary to run the
models in chunks rather than all in one go. To do this, you can break
the transcript_list_female list into parts and feed it into
the map_dfr function. The completed subset can then be
saved to your hard disk, removed from R and cleared from your computers
memory. This frees up memory to run another chunk without losing
progress. Expand the code chunk below to see an example. All other
chunks have been run and saved for later use.
# run the RFF function
transcript_list_1 <- transcript_list[1:2000]
if(!file.exists("data/transcriptome_output/transcriptome_chunks/RFF_transcript_1.csv")){
RFF_transcript_1 <- map_dfr(transcript_list_1, RFF_transcriptome_calculator)
write_csv(RFF_transcript_1,
file = "data/transcriptome_output/transcriptome_chunks/RFF_transcript_1.csv")
rm(RFF_transcript_1)
invisible(gc())
} else RFF_transcript_1 <- read_csv("data/transcriptome_output/transcriptome_chunks/RFF_transcript_1.csv")
Load all the posterior draws, combine and summarise and save the result to the hard disk. This allows us to simply load the summarised results once everything has been run once.
if(!file.exists("data/transcriptome_output/R_summarised_transcriptome.csv")){
RFF_transcriptome_complete <-
rbind(
read_csv("data/transcriptome_output/transcriptome_chunks/RFF_transcript_1.csv"),
read_csv("data/transcriptome_output/transcriptome_chunks/RFF_transcript_2.csv"),
read_csv("data/transcriptome_output/transcriptome_chunks/RFF_transcript_3.csv"),
read_csv("data/transcriptome_output/transcriptome_chunks/RFF_transcript_4.csv"),
read_csv("data/transcriptome_output/transcriptome_chunks/RFF_transcript_5.csv"),
read_csv("data/transcriptome_output/transcriptome_chunks/RFF_transcript_6.csv"),
read_csv("data/transcriptome_output/transcriptome_chunks/RFF_transcript_7.csv")
) %>%
group_by(Trait) %>%
summarise_draws("median", "sd", ~quantile(.x, probs = c(0.025, 0.975), na.rm = TRUE)) %>%
ungroup() %>%
select(-variable) %>%
mutate(absolute_R = abs(median),
Fitness_Sex = "Female",
Trait_Sex = "Female")
invisible(gc())
RMF_transcriptome_complete <-
rbind(
read_csv("data/transcriptome_output/transcriptome_chunks/RMF_transcript_1.csv"),
read_csv("data/transcriptome_output/transcriptome_chunks/RMF_transcript_2.csv"),
read_csv("data/transcriptome_output/transcriptome_chunks/RMF_transcript_3.csv"),
read_csv("data/transcriptome_output/transcriptome_chunks/RMF_transcript_4.csv"),
read_csv("data/transcriptome_output/transcriptome_chunks/RMF_transcript_5.csv"),
read_csv("data/transcriptome_output/transcriptome_chunks/RMF_transcript_6.csv"),
read_csv("data/transcriptome_output/transcriptome_chunks/RMF_transcript_7.csv")
) %>%
group_by(Trait) %>%
summarise_draws("median", "sd", ~quantile(.x, probs = c(0.025, 0.975), na.rm = TRUE)) %>%
ungroup() %>%
select(-variable) %>%
mutate(absolute_R = abs(median),
Fitness_Sex = "Male",
Trait_Sex = "Female")
invisible(gc())
RMM_transcriptome_complete <-
rbind(
read_csv("data/transcriptome_output/transcriptome_chunks/RMM_transcript_1.csv"),
read_csv("data/transcriptome_output/transcriptome_chunks/RMM_transcript_2.csv"),
read_csv("data/transcriptome_output/transcriptome_chunks/RMM_transcript_3.csv"),
read_csv("data/transcriptome_output/transcriptome_chunks/RMM_transcript_4.csv"),
read_csv("data/transcriptome_output/transcriptome_chunks/RMM_transcript_5.csv"),
read_csv("data/transcriptome_output/transcriptome_chunks/RMM_transcript_6.csv"),
read_csv("data/transcriptome_output/transcriptome_chunks/RMM_transcript_7.csv")
) %>%
group_by(Trait) %>%
summarise_draws("median", "sd", ~quantile(.x, probs = c(0.025, 0.975), na.rm = TRUE)) %>%
ungroup() %>%
select(-variable) %>%
mutate(absolute_R = abs(median),
Fitness_Sex = "Male",
Trait_Sex = "Male")
invisible(gc())
RFM_transcriptome_complete <-
rbind(
read_csv("data/transcriptome_output/transcriptome_chunks/RFM_transcript_1.csv"),
read_csv("data/transcriptome_output/transcriptome_chunks/RFM_transcript_2.csv"),
read_csv("data/transcriptome_output/transcriptome_chunks/RFM_transcript_3.csv"),
read_csv("data/transcriptome_output/transcriptome_chunks/RFM_transcript_4.csv"),
read_csv("data/transcriptome_output/transcriptome_chunks/RFM_transcript_5.csv"),
read_csv("data/transcriptome_output/transcriptome_chunks/RFM_transcript_6.csv"),
read_csv("data/transcriptome_output/transcriptome_chunks/RFM_transcript_7.csv")
) %>%
group_by(Trait) %>%
summarise_draws("median", "sd", ~quantile(.x, probs = c(0.025, 0.975), na.rm = TRUE)) %>%
ungroup() %>%
select(-variable) %>%
mutate(absolute_R = abs(median),
Fitness_Sex = "Female",
Trait_Sex = "Male")
invisible(gc())
R_summarised_transcriptome <-
bind_rows(RFF_transcriptome_complete, RMF_transcriptome_complete,
RFM_transcriptome_complete, RMM_transcriptome_complete)
write_csv(R_summarised_transcriptome, "data/transcriptome_output/R_summarised_transcriptome.csv")
} else R_summarised_transcriptome <- read_csv("data/transcriptome_output/R_summarised_transcriptome.csv")
\(~\)
First, find \(R_{overall}\)
if(!file.exists("data/transcriptome_output/R_overall_transcript_summarised.csv")){
R_overall_transcript_summarised <-
bind_rows(
left_join(read_csv(
"data/transcriptome_output/transcriptome_chunks/RFF_transcript_1.csv") %>%
group_by(Trait) %>%
mutate(draw = row_number()) %>%
ungroup(),
read_csv(
"data/transcriptome_output/transcriptome_chunks/RMF_transcript_1.csv") %>%
group_by(Trait) %>%
mutate(draw = row_number()) %>%
ungroup(), by = c("Trait", "draw")) %>%
mutate(R_overall =(Response_to_selection_female + Response_to_selection_male)/2) %>%
select(Trait, R_overall) %>%
group_by(Trait) %>%
summarise_draws("median", "sd", ~quantile(.x, probs = c(0.025, 0.975), na.rm = TRUE)) %>%
ungroup() %>%
rename(R_overall = median) %>%
select(-variable) %>%
mutate(Trait_Sex = "Female"),
left_join(read_csv(
"data/transcriptome_output/transcriptome_chunks/RFF_transcript_2.csv") %>%
group_by(Trait) %>%
mutate(draw = row_number()) %>%
ungroup(),
read_csv(
"data/transcriptome_output/transcriptome_chunks/RMF_transcript_2.csv") %>%
group_by(Trait) %>%
mutate(draw = row_number()) %>%
ungroup(), by = c("Trait", "draw")) %>%
mutate(R_overall =(Response_to_selection_female + Response_to_selection_male)/2) %>%
select(Trait, R_overall) %>%
group_by(Trait) %>%
summarise_draws("median", "sd", ~quantile(.x, probs = c(0.025, 0.975), na.rm = TRUE)) %>%
ungroup() %>%
rename(R_overall = median) %>%
select(-variable) %>%
mutate(Trait_Sex = "Female"),
left_join(read_csv(
"data/transcriptome_output/transcriptome_chunks/RFF_transcript_3.csv") %>%
group_by(Trait) %>%
mutate(draw = row_number()) %>%
ungroup(),
read_csv(
"data/transcriptome_output/transcriptome_chunks/RMF_transcript_3.csv") %>%
group_by(Trait) %>%
mutate(draw = row_number()) %>%
ungroup(), by = c("Trait", "draw")) %>%
mutate(R_overall =(Response_to_selection_female + Response_to_selection_male)/2) %>%
select(Trait, R_overall) %>%
group_by(Trait) %>%
summarise_draws("median", "sd", ~quantile(.x, probs = c(0.025, 0.975), na.rm = TRUE)) %>%
ungroup() %>%
rename(R_overall = median) %>%
select(-variable) %>%
mutate(Trait_Sex = "Female"),
left_join(read_csv(
"data/transcriptome_output/transcriptome_chunks/RFF_transcript_4.csv") %>%
group_by(Trait) %>%
mutate(draw = row_number()) %>%
ungroup(),
read_csv(
"data/transcriptome_output/transcriptome_chunks/RMF_transcript_4.csv") %>%
group_by(Trait) %>%
mutate(draw = row_number()) %>%
ungroup(), by = c("Trait", "draw")) %>%
mutate(R_overall =(Response_to_selection_female + Response_to_selection_male)/2) %>%
select(Trait, R_overall) %>%
group_by(Trait) %>%
summarise_draws("median", "sd", ~quantile(.x, probs = c(0.025, 0.975), na.rm = TRUE)) %>%
ungroup() %>%
rename(R_overall = median) %>%
select(-variable) %>%
mutate(Trait_Sex = "Female"),
left_join(read_csv(
"data/transcriptome_output/transcriptome_chunks/RFF_transcript_5.csv") %>%
group_by(Trait) %>%
mutate(draw = row_number()) %>%
ungroup(),
read_csv(
"data/transcriptome_output/transcriptome_chunks/RMF_transcript_5.csv") %>%
group_by(Trait) %>%
mutate(draw = row_number()) %>%
ungroup(), by = c("Trait", "draw")) %>%
mutate(R_overall =(Response_to_selection_female + Response_to_selection_male)/2) %>%
select(Trait, R_overall) %>%
group_by(Trait) %>%
summarise_draws("median", "sd", ~quantile(.x, probs = c(0.025, 0.975), na.rm = TRUE)) %>%
ungroup() %>%
rename(R_overall = median) %>%
select(-variable) %>%
mutate(Trait_Sex = "Female"),
left_join(read_csv(
"data/transcriptome_output/transcriptome_chunks/RFF_transcript_6.csv") %>%
group_by(Trait) %>%
mutate(draw = row_number()) %>%
ungroup(),
read_csv(
"data/transcriptome_output/transcriptome_chunks/RMF_transcript_6.csv") %>%
group_by(Trait) %>%
mutate(draw = row_number()) %>%
ungroup(), by = c("Trait", "draw")) %>%
mutate(R_overall =(Response_to_selection_female + Response_to_selection_male)/2) %>%
select(Trait, R_overall) %>%
group_by(Trait) %>%
summarise_draws("median", "sd", ~quantile(.x, probs = c(0.025, 0.975), na.rm = TRUE)) %>%
ungroup() %>%
rename(R_overall = median) %>%
select(-variable) %>%
mutate(Trait_Sex = "Female"),
left_join(read_csv(
"data/transcriptome_output/transcriptome_chunks/RFF_transcript_7.csv") %>%
group_by(Trait) %>%
mutate(draw = row_number()) %>%
ungroup(),
read_csv(
"data/transcriptome_output/transcriptome_chunks/RMF_transcript_7.csv") %>%
group_by(Trait) %>%
mutate(draw = row_number()) %>%
ungroup(), by = c("Trait", "draw")) %>%
mutate(R_overall =(Response_to_selection_female + Response_to_selection_male)/2) %>%
select(Trait, R_overall) %>%
group_by(Trait) %>%
summarise_draws("median", "sd", ~quantile(.x, probs = c(0.025, 0.975), na.rm = TRUE)) %>%
ungroup() %>%
rename(R_overall = median) %>%
select(-variable) %>%
mutate(Trait_Sex = "Female"),
# now the traits measured in males
left_join(read_csv(
"data/transcriptome_output/transcriptome_chunks/RMM_transcript_1.csv") %>%
group_by(Trait) %>%
mutate(draw = row_number()) %>%
ungroup(),
read_csv("data/transcriptome_output/transcriptome_chunks/RFM_transcript_1.csv") %>%
group_by(Trait) %>%
mutate(draw = row_number()) %>%
ungroup(), by = c("Trait", "draw")) %>%
mutate(R_overall =(Response_to_selection_female + Response_to_selection_male)/2) %>%
select(Trait, R_overall) %>%
group_by(Trait) %>%
summarise_draws("median", "sd", ~quantile(.x, probs = c(0.025, 0.975), na.rm = TRUE)) %>%
ungroup() %>%
rename(R_overall = median) %>%
select(-variable) %>%
mutate(Trait_Sex = "Male"),
left_join(read_csv("data/transcriptome_output/transcriptome_chunks/RMM_transcript_2.csv") %>%
group_by(Trait) %>%
mutate(draw = row_number()) %>%
ungroup(),
read_csv("data/transcriptome_output/transcriptome_chunks/RFM_transcript_2.csv") %>%
group_by(Trait) %>%
mutate(draw = row_number()) %>%
ungroup(), by = c("Trait", "draw")) %>%
mutate(R_overall =(Response_to_selection_female + Response_to_selection_male)/2) %>%
select(Trait, R_overall) %>%
group_by(Trait) %>%
summarise_draws("median", "sd", ~quantile(.x, probs = c(0.025, 0.975), na.rm = TRUE)) %>%
ungroup() %>%
rename(R_overall = median) %>%
select(-variable) %>%
mutate(Trait_Sex = "Male"),
left_join(read_csv("data/transcriptome_output/transcriptome_chunks/RMM_transcript_3.csv") %>%
group_by(Trait) %>%
mutate(draw = row_number()) %>%
ungroup(),
read_csv("data/transcriptome_output/transcriptome_chunks/RFM_transcript_3.csv") %>%
group_by(Trait) %>%
mutate(draw = row_number()) %>%
ungroup(), by = c("Trait", "draw")) %>%
mutate(R_overall =(Response_to_selection_female + Response_to_selection_male)/2) %>%
select(Trait, R_overall) %>%
group_by(Trait) %>%
summarise_draws("median", "sd", ~quantile(.x, probs = c(0.025, 0.975), na.rm = TRUE)) %>%
ungroup() %>%
rename(R_overall = median) %>%
select(-variable) %>%
mutate(Trait_Sex = "Male"),
left_join(read_csv("data/transcriptome_output/transcriptome_chunks/RMM_transcript_4.csv") %>%
group_by(Trait) %>%
mutate(draw = row_number()) %>%
ungroup(),
read_csv("data/transcriptome_output/transcriptome_chunks/RFM_transcript_4.csv") %>%
group_by(Trait) %>%
mutate(draw = row_number()) %>%
ungroup(), by = c("Trait", "draw")) %>%
mutate(R_overall =(Response_to_selection_female + Response_to_selection_male)/2) %>%
select(Trait, R_overall) %>%
group_by(Trait) %>%
summarise_draws("median", "sd", ~quantile(.x, probs = c(0.025, 0.975), na.rm = TRUE)) %>%
ungroup() %>%
rename(R_overall = median) %>%
select(-variable) %>%
mutate(Trait_Sex = "Male"),
left_join(read_csv("data/transcriptome_output/transcriptome_chunks/RMM_transcript_5.csv") %>%
group_by(Trait) %>%
mutate(draw = row_number()) %>%
ungroup(),
read_csv("data/transcriptome_output/transcriptome_chunks/RFM_transcript_5.csv") %>%
group_by(Trait) %>%
mutate(draw = row_number()) %>%
ungroup(), by = c("Trait", "draw")) %>%
mutate(R_overall =(Response_to_selection_female + Response_to_selection_male)/2) %>%
select(Trait, R_overall) %>%
group_by(Trait) %>%
summarise_draws("median", "sd", ~quantile(.x, probs = c(0.025, 0.975), na.rm = TRUE)) %>%
ungroup() %>%
rename(R_overall = median) %>%
select(-variable) %>%
mutate(Trait_Sex = "Male"),
left_join(read_csv("data/transcriptome_output/transcriptome_chunks/RMM_transcript_6.csv") %>%
group_by(Trait) %>%
mutate(draw = row_number()) %>%
ungroup(),
read_csv("data/transcriptome_output/transcriptome_chunks/RFM_transcript_6.csv") %>%
group_by(Trait) %>%
mutate(draw = row_number()) %>%
ungroup(), by = c("Trait", "draw")) %>%
mutate(R_overall =(Response_to_selection_female + Response_to_selection_male)/2) %>%
select(Trait, R_overall) %>%
group_by(Trait) %>%
summarise_draws("median", "sd", ~quantile(.x, probs = c(0.025, 0.975), na.rm = TRUE)) %>%
ungroup() %>%
rename(R_overall = median) %>%
select(-variable) %>%
mutate(Trait_Sex = "Male"),
left_join(read_csv("data/transcriptome_output/transcriptome_chunks/RMM_transcript_7.csv") %>%
group_by(Trait) %>%
mutate(draw = row_number()) %>%
ungroup(),
read_csv("data/transcriptome_output/transcriptome_chunks/RFM_transcript_7.csv") %>%
group_by(Trait) %>%
mutate(draw = row_number()) %>%
ungroup(), by = c("Trait", "draw")) %>%
mutate(R_overall =(Response_to_selection_female + Response_to_selection_male)/2) %>%
select(Trait, R_overall) %>%
group_by(Trait) %>%
summarise_draws("median", "sd", ~quantile(.x, probs = c(0.025, 0.975), na.rm = TRUE)) %>%
ungroup() %>%
rename(R_overall = median) %>%
select(-variable) %>%
mutate(Trait_Sex = "Male")
)
write_csv(R_overall_transcript_summarised, "data/transcriptome_output/R_overall_transcript_summarised.csv")
} else R_overall_transcript_summarised <- read_csv("data/transcriptome_output/R_overall_transcript_summarised.csv")
Now find \(\Delta R\)
if(!file.exists("data/transcriptome_output/R_diff_transcript_summarised.csv")){
R_diff_transcript_summarised <-
bind_rows(
left_join(read_csv("data/transcriptome_output/transcriptome_chunks/RFF_transcript_1.csv") %>%
group_by(Trait) %>%
mutate(draw = row_number()) %>%
ungroup(),
read_csv("data/transcriptome_output/transcriptome_chunks/RMF_transcript_1.csv") %>%
group_by(Trait) %>%
mutate(draw = row_number()) %>%
ungroup(), by = c("Trait", "draw")) %>%
mutate(R_diff = Response_to_selection_female - Response_to_selection_male) %>%
select(Trait, R_diff) %>%
group_by(Trait) %>%
summarise_draws("median", "sd", ~quantile(.x, probs = c(0.025, 0.975), na.rm = TRUE)) %>%
ungroup() %>%
rename(R_diff = median) %>%
select(-variable) %>%
mutate(Trait_Sex = "Female"),
left_join(read_csv("data/transcriptome_output/transcriptome_chunks/RFF_transcript_2.csv") %>%
group_by(Trait) %>%
mutate(draw = row_number()) %>%
ungroup(),
read_csv("data/transcriptome_output/transcriptome_chunks/RMF_transcript_2.csv") %>%
group_by(Trait) %>%
mutate(draw = row_number()) %>%
ungroup(), by = c("Trait", "draw")) %>%
mutate(R_diff = Response_to_selection_female - Response_to_selection_male) %>%
select(Trait, R_diff) %>%
group_by(Trait) %>%
summarise_draws("median", "sd", ~quantile(.x, probs = c(0.025, 0.975), na.rm = TRUE)) %>%
ungroup() %>%
rename(R_diff = median) %>%
select(-variable) %>%
mutate(Trait_Sex = "Female"),
left_join(read_csv("data/transcriptome_output/transcriptome_chunks/RFF_transcript_3.csv") %>%
group_by(Trait) %>%
mutate(draw = row_number()) %>%
ungroup(),
read_csv("data/transcriptome_output/transcriptome_chunks/RMF_transcript_3.csv") %>%
group_by(Trait) %>%
mutate(draw = row_number()) %>%
ungroup(), by = c("Trait", "draw")) %>%
mutate(R_diff = Response_to_selection_female - Response_to_selection_male) %>%
select(Trait, R_diff) %>%
group_by(Trait) %>%
summarise_draws("median", "sd", ~quantile(.x, probs = c(0.025, 0.975), na.rm = TRUE)) %>%
ungroup() %>%
rename(R_diff = median) %>%
select(-variable) %>%
mutate(Trait_Sex = "Female"),
left_join(read_csv("data/transcriptome_output/transcriptome_chunks/RFF_transcript_4.csv") %>%
group_by(Trait) %>%
mutate(draw = row_number()) %>%
ungroup(),
read_csv("data/transcriptome_output/transcriptome_chunks/RMF_transcript_4.csv") %>%
group_by(Trait) %>%
mutate(draw = row_number()) %>%
ungroup(), by = c("Trait", "draw")) %>%
mutate(R_diff = Response_to_selection_female - Response_to_selection_male) %>%
select(Trait, R_diff) %>%
group_by(Trait) %>%
summarise_draws("median", "sd", ~quantile(.x, probs = c(0.025, 0.975), na.rm = TRUE)) %>%
ungroup() %>%
rename(R_diff = median) %>%
select(-variable) %>%
mutate(Trait_Sex = "Female"),
left_join(read_csv("data/transcriptome_output/transcriptome_chunks/RFF_transcript_5.csv") %>%
group_by(Trait) %>%
mutate(draw = row_number()) %>%
ungroup(),
read_csv("data/transcriptome_output/transcriptome_chunks/RMF_transcript_5.csv") %>%
group_by(Trait) %>%
mutate(draw = row_number()) %>%
ungroup(), by = c("Trait", "draw")) %>%
mutate(R_diff = Response_to_selection_female - Response_to_selection_male) %>%
select(Trait, R_diff) %>%
group_by(Trait) %>%
summarise_draws("median", "sd", ~quantile(.x, probs = c(0.025, 0.975), na.rm = TRUE)) %>%
ungroup() %>%
rename(R_diff = median) %>%
select(-variable) %>%
mutate(Trait_Sex = "Female"),
left_join(read_csv("data/transcriptome_output/transcriptome_chunks/RFF_transcript_6.csv") %>%
group_by(Trait) %>%
mutate(draw = row_number()) %>%
ungroup(),
read_csv("data/transcriptome_output/transcriptome_chunks/RMF_transcript_6.csv") %>%
group_by(Trait) %>%
mutate(draw = row_number()) %>%
ungroup(), by = c("Trait", "draw")) %>%
mutate(R_diff = Response_to_selection_female - Response_to_selection_male) %>%
select(Trait, R_diff) %>%
group_by(Trait) %>%
summarise_draws("median", "sd", ~quantile(.x, probs = c(0.025, 0.975), na.rm = TRUE)) %>%
ungroup() %>%
rename(R_diff = median) %>%
select(-variable) %>%
mutate(Trait_Sex = "Female"),
left_join(read_csv("data/transcriptome_output/transcriptome_chunks/RFF_transcript_7.csv") %>%
group_by(Trait) %>%
mutate(draw = row_number()) %>%
ungroup(),
read_csv("data/transcriptome_output/transcriptome_chunks/RMF_transcript_7.csv") %>%
group_by(Trait) %>%
mutate(draw = row_number()) %>%
ungroup(), by = c("Trait", "draw")) %>%
mutate(R_diff = Response_to_selection_female - Response_to_selection_male) %>%
select(Trait, R_diff) %>%
group_by(Trait) %>%
summarise_draws("median", "sd", ~quantile(.x, probs = c(0.025, 0.975), na.rm = TRUE)) %>%
ungroup() %>%
rename(R_diff = median) %>%
select(-variable) %>%
mutate(Trait_Sex = "Female"),
# now the traits measured in males
left_join(read_csv("data/transcriptome_output/transcriptome_chunks/RMM_transcript_1.csv") %>%
group_by(Trait) %>%
mutate(draw = row_number()) %>%
ungroup(),
read_csv("data/transcriptome_output/transcriptome_chunks/RFM_transcript_1.csv") %>%
group_by(Trait) %>%
mutate(draw = row_number()) %>%
ungroup(), by = c("Trait", "draw")) %>%
mutate(R_diff = Response_to_selection_female - Response_to_selection_male) %>%
select(Trait, R_diff) %>%
group_by(Trait) %>%
summarise_draws("median", "sd", ~quantile(.x, probs = c(0.025, 0.975), na.rm = TRUE)) %>%
ungroup() %>%
rename(R_diff = median) %>%
select(-variable) %>%
mutate(Trait_Sex = "Male"),
left_join(read_csv("data/transcriptome_output/transcriptome_chunks/RMM_transcript_2.csv") %>%
group_by(Trait) %>%
mutate(draw = row_number()) %>%
ungroup(),
read_csv("data/transcriptome_output/transcriptome_chunks/RFM_transcript_2.csv") %>%
group_by(Trait) %>%
mutate(draw = row_number()) %>%
ungroup(), by = c("Trait", "draw")) %>%
mutate(R_diff = Response_to_selection_female - Response_to_selection_male) %>%
select(Trait, R_diff) %>%
group_by(Trait) %>%
summarise_draws("median", "sd", ~quantile(.x, probs = c(0.025, 0.975), na.rm = TRUE)) %>%
ungroup() %>%
rename(R_diff = median) %>%
select(-variable) %>%
mutate(Trait_Sex = "Male"),
left_join(read_csv("data/transcriptome_output/transcriptome_chunks/RMM_transcript_3.csv") %>%
group_by(Trait) %>%
mutate(draw = row_number()) %>%
ungroup(),
read_csv("data/transcriptome_output/transcriptome_chunks/RFM_transcript_3.csv") %>%
group_by(Trait) %>%
mutate(draw = row_number()) %>%
ungroup(), by = c("Trait", "draw")) %>%
mutate(R_diff = Response_to_selection_female - Response_to_selection_male) %>%
select(Trait, R_diff) %>%
group_by(Trait) %>%
summarise_draws("median", "sd", ~quantile(.x, probs = c(0.025, 0.975), na.rm = TRUE)) %>%
ungroup() %>%
rename(R_diff = median) %>%
select(-variable) %>%
mutate(Trait_Sex = "Male"),
left_join(read_csv("data/transcriptome_output/transcriptome_chunks/RMM_transcript_4.csv") %>%
group_by(Trait) %>%
mutate(draw = row_number()) %>%
ungroup(),
read_csv("data/transcriptome_output/transcriptome_chunks/RFM_transcript_4.csv") %>%
group_by(Trait) %>%
mutate(draw = row_number()) %>%
ungroup(), by = c("Trait", "draw")) %>%
mutate(R_diff = Response_to_selection_female - Response_to_selection_male) %>%
select(Trait, R_diff) %>%
group_by(Trait) %>%
summarise_draws("median", "sd", ~quantile(.x, probs = c(0.025, 0.975), na.rm = TRUE)) %>%
ungroup() %>%
rename(R_diff = median) %>%
select(-variable) %>%
mutate(Trait_Sex = "Male"),
left_join(read_csv("data/transcriptome_output/transcriptome_chunks/RMM_transcript_5.csv") %>%
group_by(Trait) %>%
mutate(draw = row_number()) %>%
ungroup(),
read_csv("data/transcriptome_output/transcriptome_chunks/RFM_transcript_5.csv") %>%
group_by(Trait) %>%
mutate(draw = row_number()) %>%
ungroup(), by = c("Trait", "draw")) %>%
mutate(R_diff = Response_to_selection_female - Response_to_selection_male) %>%
select(Trait, R_diff) %>%
group_by(Trait) %>%
summarise_draws("median", "sd", ~quantile(.x, probs = c(0.025, 0.975), na.rm = TRUE)) %>%
ungroup() %>%
rename(R_diff = median) %>%
select(-variable) %>%
mutate(Trait_Sex = "Male"),
left_join(read_csv("data/transcriptome_output/transcriptome_chunks/RMM_transcript_6.csv") %>%
group_by(Trait) %>%
mutate(draw = row_number()) %>%
ungroup(),
read_csv("data/transcriptome_output/transcriptome_chunks/RFM_transcript_6.csv") %>%
group_by(Trait) %>%
mutate(draw = row_number()) %>%
ungroup(), by = c("Trait", "draw")) %>%
mutate(R_diff = Response_to_selection_female - Response_to_selection_male) %>%
select(Trait, R_diff) %>%
group_by(Trait) %>%
summarise_draws("median", "sd", ~quantile(.x, probs = c(0.025, 0.975), na.rm = TRUE)) %>%
ungroup() %>%
rename(R_diff = median) %>%
select(-variable) %>%
mutate(Trait_Sex = "Male"),
left_join(read_csv("data/transcriptome_output/transcriptome_chunks/RMM_transcript_7.csv") %>%
group_by(Trait) %>%
mutate(draw = row_number()) %>%
ungroup(),
read_csv("data/transcriptome_output/transcriptome_chunks/RFM_transcript_7.csv") %>%
group_by(Trait) %>%
mutate(draw = row_number()) %>%
ungroup(), by = c("Trait", "draw")) %>%
mutate(R_diff = Response_to_selection_female - Response_to_selection_male) %>%
select(Trait, R_diff) %>%
group_by(Trait) %>%
summarise_draws("median", "sd", ~quantile(.x, probs = c(0.025, 0.975), na.rm = TRUE)) %>%
ungroup() %>%
rename(R_diff = median) %>%
select(-variable) %>%
mutate(Trait_Sex = "Male")
)
write_csv(R_diff_transcript_summarised, "data/transcriptome_output/R_diff_transcript_summarised.csv")
} else R_diff_transcript_summarised <- read_csv("data/transcriptome_output/R_diff_transcript_summarised.csv")
Use the brms hypothesis() function to
compute evidence ratios, just as we did for the organismal level
phenotypic traits.
# build a function to calculate the evidence ratio, get the relevant info from the output and convert it to a tibble. Basically the same function, just with a different dataframe input and a memory clearing step that slows it down but allows it to run for many traits.
find_evidence_ratios_transcriptome <- function(Trait_name, specify_hypothesis){
x <- hypothesis(R_transcriptome %>%
filter(Trait == Trait_name),
specify_hypothesis)
x <- x$hypothesis
invisible(gc())
x %>% as_tibble() %>%
select(Evid.Ratio) %>%
mutate(Trait = Trait_name)
}
This pc has 16gb memory; not enough to load in all the transcript
posterior draws at once. Having previously saved these to the hard disk
we load the draws for all female measured traits and calculate the
evidence ratio that \(R_{FF} > 0\),
\(R_{MF} > 0\), \(R_F > 0\) and \(\Delta R > 0\). The draws for the female
measured traits were then removed from the computers memory using the
rm() and gc() functions and the process was
repeated for the male measured traits.
# find evidence ratios for the traits measured in females
if(!file.exists("data/transcriptome_output/evidence_ratios_transcriptome_female.csv")){ # name refers to traits measured in females
R_transcriptome <-
cbind(
rbind(
read_csv("data/transcriptome_output/transcriptome_chunks/RFF_transcript_1.csv"),
read_csv("data/transcriptome_output/transcriptome_chunks/RFF_transcript_2.csv"),
read_csv("data/transcriptome_output/transcriptome_chunks/RFF_transcript_3.csv"),
read_csv("data/transcriptome_output/transcriptome_chunks/RFF_transcript_4.csv"),
read_csv("data/transcriptome_output/transcriptome_chunks/RFF_transcript_5.csv"),
read_csv("data/transcriptome_output/transcriptome_chunks/RFF_transcript_6.csv"),
read_csv("data/transcriptome_output/transcriptome_chunks/RFF_transcript_7.csv")
),
rbind(
read_csv("data/transcriptome_output/transcriptome_chunks/RMF_transcript_1.csv"),
read_csv("data/transcriptome_output/transcriptome_chunks/RMF_transcript_2.csv"),
read_csv("data/transcriptome_output/transcriptome_chunks/RMF_transcript_3.csv"),
read_csv("data/transcriptome_output/transcriptome_chunks/RMF_transcript_4.csv"),
read_csv("data/transcriptome_output/transcriptome_chunks/RMF_transcript_5.csv"),
read_csv("data/transcriptome_output/transcriptome_chunks/RMF_transcript_6.csv"),
read_csv("data/transcriptome_output/transcriptome_chunks/RMF_transcript_7.csv")) %>%
select(-Trait)
)
R_transcriptome <-
R_transcriptome %>%
mutate(R_overall = (Response_to_selection_female - Response_to_selection_male) / 2,
R_diff = Response_to_selection_female - Response_to_selection_male)
# calculate evidence ratios
R_evidence_ratios_female <-
map_dfr(transcript_list,
find_evidence_ratios_transcriptome,
"Response_to_selection_female > 0") %>%
rename(Female_Evid_Ratio = Evid.Ratio)
R_evidence_ratios_male <-
map_dfr(transcript_list,
find_evidence_ratios_transcriptome,
"Response_to_selection_male > 0") %>%
rename(Male_Evid_Ratio = Evid.Ratio)
R_evidence_ratios_overall <-
map_dfr(transcript_list,
find_evidence_ratios_transcriptome,
"R_overall > 0") %>%
rename(Overall_Evid_Ratio = Evid.Ratio)
R_evidence_ratios_diff <-
map_dfr(transcript_list,
find_evidence_ratios_transcriptome,
"R_diff > 0") %>%
rename(Diff_Evid_Ratio = Evid.Ratio)
# write each result to csv in
evidence_ratios_transcriptome_female <-
R_evidence_ratios_female %>%
left_join(R_evidence_ratios_male) %>%
left_join(R_evidence_ratios_overall) %>%
left_join(R_evidence_ratios_diff) %>%
select(Trait, everything()) %>%
rename(Female_Evid_Ratio = Female_evidence_ratio) %>% # fixing small typo
mutate(Trait_Sex = "Female")
write_csv(evidence_ratios_transcriptome_female, "data/transcriptome_output/evidence_ratios_transcriptome_female.csv")
} else evidence_ratios_transcriptome_female <- read_csv("data/transcriptome_output/evidence_ratios_transcriptome_female.csv")
# find evidence ratios for the traits measured in males
if(!file.exists("data/transcriptome_output/evidence_ratios_transcriptome_male.csv")){ # name refers to traits measured in females
R_transcriptome <-
cbind(
rbind(
read_csv("data/transcriptome_output/transcriptome_chunks/RFM_transcript_1.csv"),
read_csv("data/transcriptome_output/transcriptome_chunks/RFM_transcript_2.csv"),
read_csv("data/transcriptome_output/transcriptome_chunks/RFM_transcript_3.csv"),
read_csv("data/transcriptome_output/transcriptome_chunks/RFM_transcript_4.csv"),
read_csv("data/transcriptome_output/transcriptome_chunks/RFM_transcript_5.csv"),
read_csv("data/transcriptome_output/transcriptome_chunks/RFM_transcript_6.csv"),
read_csv("data/transcriptome_output/transcriptome_chunks/RFM_transcript_7.csv")
),
rbind(
read_csv("data/transcriptome_output/transcriptome_chunks/RMM_transcript_1.csv"),
read_csv("data/transcriptome_output/transcriptome_chunks/RMM_transcript_2.csv"),
read_csv("data/transcriptome_output/transcriptome_chunks/RMM_transcript_3.csv"),
read_csv("data/transcriptome_output/transcriptome_chunks/RMM_transcript_4.csv"),
read_csv("data/transcriptome_output/transcriptome_chunks/RMM_transcript_5.csv"),
read_csv("data/transcriptome_output/transcriptome_chunks/RMM_transcript_6.csv"),
read_csv("data/transcriptome_output/transcriptome_chunks/RMM_transcript_7.csv")) %>%
select(-Trait)
)
R_transcriptome <-
R_transcriptome %>%
mutate(R_overall = (Response_to_selection_female - Response_to_selection_male) / 2,
R_diff = Response_to_selection_female - Response_to_selection_male)
# calculate evidence ratios
R_evidence_ratios_female <-
map_dfr(transcript_list,
find_evidence_ratios_transcriptome,
"Response_to_selection_female > 0") %>%
rename(Female_Evid_Ratio = Evid.Ratio)
R_evidence_ratios_male <-
map_dfr(transcript_list,
find_evidence_ratios_transcriptome,
"Response_to_selection_male > 0") %>%
rename(Male_Evid_Ratio = Evid.Ratio)
R_evidence_ratios_overall <-
map_dfr(transcript_list,
find_evidence_ratios_transcriptome,
"R_overall > 0") %>%
rename(Overall_Evid_Ratio = Evid.Ratio)
R_evidence_ratios_diff <-
map_dfr(transcript_list,
find_evidence_ratios_transcriptome,
"R_diff > 0") %>%
rename(Diff_Evid_Ratio = Evid.Ratio)
evidence_ratios_transcriptome_male <-
R_evidence_ratios_female %>%
left_join(R_evidence_ratios_male) %>%
left_join(R_evidence_ratios_overall) %>%
left_join(R_evidence_ratios_diff) %>%
select(Trait, everything()) %>%
rename(Female_Evid_Ratio = Female_evidence_ratio) %>% # fixing small typo
mutate(Trait_Sex = "Male")
write_csv(evidence_ratios_transcriptome_male, "data/transcriptome_output/evidence_ratios_transcriptome_male.csv")
} else evidence_ratios_transcriptome_male <- read_csv("data/transcriptome_output/evidence_ratios_transcriptome_male.csv")
To calculate evidence ratios for sexually concordant responses to selection we use the same function as for the organismal-level phenotypic traits. Run the function in chunks, save them to hard-disk, clear memory and repeat. If already done, just load.
if(!file.exists("data/transcriptome_output/concord_evidence_ratio_data.csv")){
female_transcripts_1 <-
left_join(read_csv("data/transcriptome_output/transcriptome_chunks/RFF_transcript_1.csv") %>%
group_by(Trait) %>%
mutate(draw = row_number()) %>%
ungroup(),
read_csv("data/transcriptome_output/transcriptome_chunks/RMF_transcript_1.csv") %>%
group_by(Trait) %>%
mutate(draw = row_number()) %>%
ungroup(), by = c("Trait", "draw"))
trait_list_1 <- unique(female_transcripts_1$Trait)
f1 <- get_evidence_ratios(female_transcripts_1, trait_list_1, "Female")
rm(female_transcripts_1)
invisible(gc())
female_transcripts_2 <-
left_join(read_csv("data/transcriptome_output/transcriptome_chunks/RFF_transcript_2.csv") %>%
group_by(Trait) %>%
mutate(draw = row_number()) %>%
ungroup(),
read_csv("data/transcriptome_output/transcriptome_chunks/RMF_transcript_2.csv") %>%
group_by(Trait) %>%
mutate(draw = row_number()) %>%
ungroup(), by = c("Trait", "draw"))
trait_list_2 <- unique(female_transcripts_2$Trait)
f2 <- get_evidence_ratios(female_transcripts_2, trait_list_2, "Female")
rm(female_transcripts_2)
invisible(gc())
female_transcripts_3 <-
left_join(read_csv("data/transcriptome_output/transcriptome_chunks/RFF_transcript_3.csv") %>%
group_by(Trait) %>%
mutate(draw = row_number()) %>%
ungroup(),
read_csv("data/transcriptome_output/transcriptome_chunks/RMF_transcript_3.csv") %>%
group_by(Trait) %>%
mutate(draw = row_number()) %>%
ungroup(), by = c("Trait", "draw"))
trait_list_3 <- unique(female_transcripts_3$Trait)
f3 <- get_evidence_ratios(female_transcripts_3, trait_list_3, "Female")
rm(female_transcripts_3)
invisible(gc())
female_transcripts_4 <-
left_join(read_csv("data/transcriptome_output/transcriptome_chunks/RFF_transcript_4.csv") %>%
group_by(Trait) %>%
mutate(draw = row_number()) %>%
ungroup(),
read_csv("data/transcriptome_output/transcriptome_chunks/RMF_transcript_4.csv") %>%
group_by(Trait) %>%
mutate(draw = row_number()) %>%
ungroup(), by = c("Trait", "draw"))
trait_list_4 <- unique(female_transcripts_4$Trait)
f4 <- get_evidence_ratios(female_transcripts_4, trait_list_4, "Female")
rm(female_transcripts_4)
invisible(gc())
female_transcripts_5 <-
left_join(read_csv("data/transcriptome_output/transcriptome_chunks/RFF_transcript_5.csv") %>%
group_by(Trait) %>%
mutate(draw = row_number()) %>%
ungroup(),
read_csv("data/transcriptome_output/transcriptome_chunks/RMF_transcript_5.csv") %>%
group_by(Trait) %>%
mutate(draw = row_number()) %>%
ungroup(), by = c("Trait", "draw"))
trait_list_5 <- unique(female_transcripts_5$Trait)
f5 <- get_evidence_ratios(female_transcripts_5, trait_list_5, "Female")
rm(female_transcripts_5)
invisible(gc())
female_transcripts_6 <-
left_join(read_csv("data/transcriptome_output/transcriptome_chunks/RFF_transcript_6.csv") %>%
group_by(Trait) %>%
mutate(draw = row_number()) %>%
ungroup(),
read_csv("data/transcriptome_output/transcriptome_chunks/RMF_transcript_6.csv") %>%
group_by(Trait) %>%
mutate(draw = row_number()) %>%
ungroup(), by = c("Trait", "draw"))
trait_list_6 <- unique(female_transcripts_6$Trait)
f6 <- get_evidence_ratios(female_transcripts_6, trait_list_6, "Female")
rm(female_transcripts_6)
invisible(gc())
female_transcripts_7 <-
left_join(read_csv("data/transcriptome_output/transcriptome_chunks/RFF_transcript_7.csv") %>%
group_by(Trait) %>%
mutate(draw = row_number()) %>%
ungroup(),
read_csv("data/transcriptome_output/transcriptome_chunks/RMF_transcript_7.csv") %>%
group_by(Trait) %>%
mutate(draw = row_number()) %>%
ungroup(), by = c("Trait", "draw"))
trait_list_7 <- unique(female_transcripts_7$Trait)
f7 <- get_evidence_ratios(female_transcripts_7, trait_list_7, "Female")
rm(female_transcripts_7)
invisible(gc())
# now the transcripts measured in males
male_transcripts_1 <-
left_join(read_csv("data/transcriptome_output/transcriptome_chunks/RMM_transcript_1.csv") %>%
group_by(Trait) %>%
mutate(draw = row_number()) %>%
ungroup(),
read_csv("data/transcriptome_output/transcriptome_chunks/RFM_transcript_1.csv") %>%
group_by(Trait) %>%
mutate(draw = row_number()) %>%
ungroup(), by = c("Trait", "draw"))
m_trait_list_1 <- unique(male_transcripts_1$Trait)
m1 <- get_evidence_ratios(male_transcripts_1, m_trait_list_1, "Male")
rm(male_transcripts_1)
invisible(gc())
male_transcripts_2 <-
left_join(read_csv("data/transcriptome_output/transcriptome_chunks/RMM_transcript_2.csv") %>%
group_by(Trait) %>%
mutate(draw = row_number()) %>%
ungroup(),
read_csv("data/transcriptome_output/transcriptome_chunks/RFM_transcript_2.csv") %>%
group_by(Trait) %>%
mutate(draw = row_number()) %>%
ungroup(), by = c("Trait", "draw"))
m_trait_list_2 <- unique(male_transcripts_2$Trait)
m2 <- get_evidence_ratios(male_transcripts_2, m_trait_list_2, "Male")
rm(male_transcripts_2)
invisible(gc())
male_transcripts_3 <-
left_join(read_csv("data/transcriptome_output/transcriptome_chunks/RMM_transcript_3.csv") %>%
group_by(Trait) %>%
mutate(draw = row_number()) %>%
ungroup(),
read_csv("data/transcriptome_output/transcriptome_chunks/RFM_transcript_3.csv") %>%
group_by(Trait) %>%
mutate(draw = row_number()) %>%
ungroup(), by = c("Trait", "draw"))
m_trait_list_3 <- unique(male_transcripts_3$Trait)
m3 <- get_evidence_ratios(male_transcripts_3, m_trait_list_3, "Male")
rm(male_transcripts_3)
invisible(gc())
male_transcripts_4 <-
left_join(read_csv("data/transcriptome_output/transcriptome_chunks/RMM_transcript_4.csv") %>%
group_by(Trait) %>%
mutate(draw = row_number()) %>%
ungroup(),
read_csv("data/transcriptome_output/transcriptome_chunks/RFM_transcript_4.csv") %>%
group_by(Trait) %>%
mutate(draw = row_number()) %>%
ungroup(), by = c("Trait", "draw"))
m_trait_list_4 <- unique(male_transcripts_4$Trait)
m4 <- get_evidence_ratios(male_transcripts_4, m_trait_list_4, "Male")
rm(male_transcripts_4)
invisible(gc())
male_transcripts_5 <-
left_join(read_csv("data/transcriptome_output/transcriptome_chunks/RMM_transcript_5.csv") %>%
group_by(Trait) %>%
mutate(draw = row_number()) %>%
ungroup(),
read_csv("data/transcriptome_output/transcriptome_chunks/RFM_transcript_5.csv") %>%
group_by(Trait) %>%
mutate(draw = row_number()) %>%
ungroup(), by = c("Trait", "draw"))
m_trait_list_5 <- unique(male_transcripts_5$Trait)
m5 <- get_evidence_ratios(male_transcripts_5, m_trait_list_5, "Male")
rm(male_transcripts_5)
invisible(gc())
male_transcripts_6 <-
left_join(read_csv("data/transcriptome_output/transcriptome_chunks/RMM_transcript_6.csv") %>%
group_by(Trait) %>%
mutate(draw = row_number()) %>%
ungroup(),
read_csv("data/transcriptome_output/transcriptome_chunks/RFM_transcript_6.csv") %>%
group_by(Trait) %>%
mutate(draw = row_number()) %>%
ungroup(), by = c("Trait", "draw"))
m_trait_list_6 <- unique(male_transcripts_6$Trait)
m6 <- get_evidence_ratios(male_transcripts_6, m_trait_list_6, "Male")
rm(male_transcripts_6)
invisible(gc())
male_transcripts_7 <-
left_join(read_csv("data/transcriptome_output/transcriptome_chunks/RMM_transcript_7.csv") %>%
group_by(Trait) %>%
mutate(draw = row_number()) %>%
ungroup(),
read_csv("data/transcriptome_output/transcriptome_chunks/RFM_transcript_7.csv") %>%
group_by(Trait) %>%
mutate(draw = row_number()) %>%
ungroup(), by = c("Trait", "draw"))
m_trait_list_7 <- unique(male_transcripts_7$Trait)
m7 <- get_evidence_ratios(male_transcripts_7, m_trait_list_7, "Male")
rm(male_transcripts_7)
invisible(gc())
concord_evidence_ratio_data <-
bind_rows(f1, f2, f3, f4, f5, f6, f7, m1, m2, m3, m4, m5, m6, m7)
write_csv(concord_evidence_ratio_data, "data/transcriptome_output/concord_evidence_ratio_data.csv")
} else concord_evidence_ratio_data <- read_csv("data/transcriptome_output/concord_evidence_ratio_data.csv")
Bind all the evidence ratios together
evidence_ratios_transcriptome <-
evidence_ratios_transcriptome_female %>%
bind_rows(evidence_ratios_transcriptome_male %>% rename(Male_Evid_Ratio = Male_evidence_ratio)) %>% # fix a typo
left_join(concord_evidence_ratio_data %>% rename(Concord_Evid_Ratio = evidence_ratio_concord) %>% select(Trait, Trait_Sex, Concord_Evid_Ratio),
by = c("Trait", "Trait_Sex")) %>%
left_join(gene_annotations %>% rename(Trait = FBID)) %>%
mutate(chromosome_numeric = case_when(chromosome == "2L" ~ 1,
chromosome == "2R" ~ 2,
chromosome == "3L" ~ 3,
chromosome == "3R" ~ 4,
chromosome == "X" ~ 5)) %>%
#left_join(R_summarised_transcriptome %>% distinct(Trait, Trait_Sex)) %>%
#mutate(Trait = fct_reorder(Trait, chromosome_numeric))
arrange(chromosome, Trait) %>%
mutate(position = 1:n(),
Trait_Sex = if_else(Trait_Sex == "Female", "Traits measured in females",
"Traits measured in males")) %>%
select(Trait, everything())
\(~\)
Find the number of evidence ratios > 39 or < 1/39 for each partition of \(R\)
response_pos <- evidence_ratios_transcriptome %>% filter(Overall_Evid_Ratio > 39) %>% arrange(-Overall_Evid_Ratio)
response_neg <- evidence_ratios_transcriptome %>% filter(Overall_Evid_Ratio < 1/39) %>% arrange(Overall_Evid_Ratio)
response_to_female_pos <- evidence_ratios_transcriptome %>% filter(Female_Evid_Ratio > 39) %>% arrange(-Female_Evid_Ratio)
response_to_female_neg <- evidence_ratios_transcriptome %>% filter(Female_Evid_Ratio < 1/39) %>% arrange(Female_Evid_Ratio)
response_to_male_pos <- evidence_ratios_transcriptome %>% filter(Male_Evid_Ratio > 39) %>% arrange(-Male_Evid_Ratio)
response_to_male_neg <- evidence_ratios_transcriptome %>% filter(Male_Evid_Ratio < 1/39) %>% arrange(Male_Evid_Ratio)
response_sexes_diff <-
evidence_ratios_transcriptome %>% filter(Diff_Evid_Ratio > 39 | Diff_Evid_Ratio < 1/39) %>% arrange(-Diff_Evid_Ratio)
sexually_concordant <- evidence_ratios_transcriptome %>% filter(Concord_Evid_Ratio > 39) %>% arrange(-Concord_Evid_Ratio)
sexually_antagonistic <- evidence_ratios_transcriptome %>% filter(Concord_Evid_Ratio < 1/39) %>% arrange(Concord_Evid_Ratio)
684 transcripts are predicted to have a positive overall response to selection, while 844 transcripts ar expected to have a negative overall response.
In response to selection on females, 424 transcripts are expected to increase in expression, whereas 530 are expected to decrease in mean expression.
In response to selection on males, 468 transcripts are expected to increase in expression, whereas 619 are expected to decrease in mean expression.
575 transcripts have notably different responses to selection on females compared with males.
1 transcripts are expected to have a sexually antagonistic response to selection, while 25 transcripts are expected to have a sexually concordant response. From these numbers, we can roughly estimate that 3.8% of transcripts that are responding to selection in both sexes are sexually antagonistic.
Build the evidence ratio plots
x_labels <- evidence_ratios_transcriptome %>%
group_by(chromosome) %>%
summarise(position = min(position) + (max(position) - min(position))/2)
guild_pal <- c("#6e7cb9", "#7bbcd5", "#d0e2af", "#f5db99",
"#e89c81", "#d2848d")
plot_transcriptome_a <-
evidence_ratios_transcriptome %>%
ggplot(aes(x = position, y = log2(Overall_Evid_Ratio))) +
geom_point(aes(colour = chromosome),
size = 2.5, alpha = 0.8, shape = 17) +
geom_hline(yintercept = log2(1), linetype = 2, linewidth = 1) +
geom_hline(yintercept = log2(39), linetype = 3, linewidth = 0.9) +
geom_hline(yintercept = log2(1/39), linetype = 3, linewidth = 0.9) +
#scale_shape_manual(values = c(25, 24)) +
scale_colour_manual(values = guild_pal) +
scale_x_continuous(breaks = x_labels$position, labels = x_labels$chromosome) +
scale_y_continuous(expand = c(0.01, 0), breaks = c(-14, -12, -10, -8, -6, -4, -2, 0, 2, 4, 6, 8, 10, 12, 14),
labels = c(paste("1/",2 ^ c(14, 12, 10, 8, 6, 4, 2), sep = ""), 2 ^ c(0, 2,4,6,8,10,12,14)), limits = c(-18, 18)) +
labs(x = NULL,
#y = "R~F~ log2(ER)",
y = "_R_ log~2~(ER)") +
geom_label_repel(data = evidence_ratios_transcriptome %>%
filter(Overall_Evid_Ratio > 1000 | Overall_Evid_Ratio < 1/1000 ),
aes(label = gene_symbol),
fill = "white",
size = 3, alpha = 0.9,
min.segment.length = 0, seed = 1,
box.padding = 0.5, point.padding = 0.5,
max.overlaps = 30) +
theme_minimal() +
theme(legend.position = "none",
panel.grid.major.x = element_blank(),
panel.grid.minor.x = element_blank(),
panel.grid.minor.y = element_blank(),
axis.text.x=element_blank(),
axis.ticks.x=element_blank(),
axis.title.x = element_text(size = 15),
axis.title.y = element_markdown(size = 10))
plot_transcriptome_b <-
evidence_ratios_transcriptome %>%
ggplot(aes(x = position, y = log2(Female_Evid_Ratio))) +
geom_point(aes(colour = chromosome),
size = 2.5, alpha = 0.8, shape = 17) +
geom_hline(yintercept = log2(1), linetype = 2, linewidth = 1) +
geom_hline(yintercept = log2(39), linetype = 3, linewidth = 0.9) +
geom_hline(yintercept = log2(1/39), linetype = 3, linewidth = 0.9) +
#scale_shape_manual(values = c(25, 24)) +
scale_colour_manual(values = guild_pal) +
scale_x_continuous(breaks = x_labels$position, labels = x_labels$chromosome) +
scale_y_continuous(expand = c(0.01, 0), breaks = c(-14, -12, -10, -8, -6, -4, -2, 0, 2, 4, 6, 8, 10, 12, 14),
labels = c(paste("1/",2 ^ c(14, 12, 10, 8, 6, 4, 2), sep = ""), 2 ^ c(0, 2,4,6,8,10,12,14)), limits = c(-16, 16)) +
labs(x = NULL,
#y = "R~F~ log2(ER)",
y = "_R_ on females log~2~(ER)") +
geom_label_repel(data = evidence_ratios_transcriptome %>%
filter(Female_Evid_Ratio > 1000 | Female_Evid_Ratio < 1/1000 ),
aes(label = gene_symbol),
fill = "white",
size = 3, alpha = 0.9,
min.segment.length = 0, seed = 1,
box.padding = 0.5, point.padding = 0.5,
max.overlaps = 30) +
theme_minimal() +
theme(legend.position = "none",
panel.grid.major.x = element_blank(),
panel.grid.minor.x = element_blank(),
panel.grid.minor.y = element_blank(),
axis.text.x=element_blank(),
axis.ticks.x=element_blank(),
axis.title.x = element_text(size = 15),
axis.title.y = element_markdown(size = 10))
plot_transcriptome_c <-
evidence_ratios_transcriptome %>%
ggplot(aes(x = position, y = log2(Male_Evid_Ratio))) +
geom_point(aes(colour = chromosome),
size = 2.5, alpha = 0.8, shape = 17) +
geom_hline(yintercept = log2(1), linetype = 2, linewidth = 1) +
geom_hline(yintercept = log2(39), linetype = 3, linewidth = 0.9) +
geom_hline(yintercept = log2(1/39), linetype = 3, linewidth = 0.9) +
#scale_shape_manual(values = c(25, 24)) +
scale_colour_manual(values = guild_pal) +
scale_x_continuous(breaks = x_labels$position, labels = x_labels$chromosome) +
scale_y_continuous(expand = c(0.01, 0), breaks = c(-14, -12, -10, -8, -6, -4, -2, 0, 2, 4, 6, 8, 10, 12, 14),
labels = c(paste("1/",2 ^ c(14, 12, 10, 8, 6, 4, 2), sep = ""), 2 ^ c(0, 2,4,6,8,10,12,14)), limits = c(-15, 15)) +
labs(x = "Chromosome arm",
#y = "R~F~ log2(ER)",
y = "_R_ on males log~2~(ER)") +
geom_label_repel(data = evidence_ratios_transcriptome %>%
filter(Male_Evid_Ratio > 1000 | Male_Evid_Ratio < 1/1000 ),
aes(label = gene_symbol),
fill = "white",
size = 3, alpha = 0.9,
min.segment.length = 0, seed = 1,
box.padding = 0.5, point.padding = 0.5,
max.overlaps = 30) +
theme_minimal() +
theme(legend.position = "none",
panel.grid.major.x = element_blank(),
panel.grid.minor.x = element_blank(),
panel.grid.minor.y = element_blank(),
axis.text.x=element_text(size = 15),
#axis.ticks.x=element_blank(),
axis.title.x = element_text(size = 15),
axis.title.y = element_markdown(size = 10))
(Figure_S4 <- plot_transcriptome_a / plot_transcriptome_b / plot_transcriptome_c +
plot_annotation(tag_levels = 'a')
)

| Version | Author | Date |
|---|---|---|
| ad196cf | tomkeaney | 2024-04-11 |
Figure S4. Triangles show evidence ratios (ERs) for gene expression traits. a shows overall \(R\), b \(R\) on females (\(R_{FF}\) and \(R_{MF}\)) and c \(R\) on males (\(R_{FM}\) and \(R_{MM}\)). The dashed lines indicate an evidence ratio of 1, where the probability that \(R\) > 0 is equal to the probability that \(R\) < 0. ER > 1 indicates a positive response is more likely, whereas ER < 1 indicates a negative response is more likely. The upper dotted line on each plot indicates an evidence ratio of 39, while the lower dotted line indicates an evidence ratio of 1/39; these correspond strongly with the frequentist \(p\) = 0.05. Traits with the most extreme evidence ratios are labelled.
plot_transcriptome_d <-
evidence_ratios_transcriptome %>%
ggplot(aes(x = position, y = log2(Diff_Evid_Ratio))) +
geom_point(aes(colour = chromosome),
size = 2.5, alpha = 0.8, shape = 17) +
geom_hline(yintercept = log2(1), linetype = 2, linewidth = 1) +
geom_hline(yintercept = log2(39), linetype = 3, linewidth = 0.9) +
geom_hline(yintercept = log2(1/39), linetype = 3, linewidth = 0.9) +
#scale_shape_manual(values = c(25, 24)) +
scale_colour_manual(values = guild_pal) +
scale_x_continuous(breaks = x_labels$position, labels = x_labels$chromosome) +
scale_y_continuous(expand = c(0.01, 0), breaks = c(-12, -10, -8, -6, -4, -2, 0, 2, 4, 6, 8, 10, 12),
labels = c(paste("1/",2 ^ c(12, 10, 8, 6, 4, 2), sep = ""), 2 ^ c(0, 2,4,6,8,10,12)), limits = c(-13, 13)) +
labs(x = "Chromosome arm",
#y = "R~F~ log2(ER)",
y = "Δ_R_ log~2~(ER)") +
geom_label_repel(data = evidence_ratios_transcriptome %>%
filter(Diff_Evid_Ratio > 250 | Diff_Evid_Ratio < 1/250 ),
aes(label = gene_symbol),
fill = "white",
size = 3, alpha = 0.9,
min.segment.length = 0, seed = 1,
box.padding = 0.5, point.padding = 0.5,
max.overlaps = 30) +
theme_minimal() +
theme(legend.position = "none",
panel.grid.major.x = element_blank(),
panel.grid.minor.x = element_blank(),
panel.grid.minor.y = element_blank(),
axis.text.x=element_blank(),
axis.ticks.x=element_blank(),
axis.title.x = element_blank(),
axis.title.y = element_markdown(size = 15))
plot_transcriptome_e <-
evidence_ratios_transcriptome %>%
ggplot(aes(x = position, y = log2(Concord_Evid_Ratio))) +
geom_point(aes(colour = chromosome),
size = 2.5, alpha = 0.8, shape = 17) +
geom_hline(yintercept = log2(1), linetype = 2, linewidth = 1) +
geom_hline(yintercept = log2(39), linetype = 3, linewidth = 0.9) +
geom_hline(yintercept = log2(1/39), linetype = 3, linewidth = 0.9) +
#scale_shape_manual(values = c(25, 24)) +
scale_colour_manual(values = guild_pal) +
scale_x_continuous(breaks = x_labels$position, labels = x_labels$chromosome) +
scale_y_continuous(expand = c(0.01, 0), breaks = c(-12, -10, -8, -6, -4, -2, 0, 2, 4, 6, 8, 10, 12),
labels = c(paste("1/",2 ^ c(12, 10, 8, 6, 4, 2), sep = ""), 2 ^ c(0, 2,4,6,8,10,12)), limits = c(-12, 12)) +
labs(x = "Chromosome arm",
y = "Sexual concrdoance log~2~(ER)") +
geom_label_repel(data = evidence_ratios_transcriptome %>%
filter(Concord_Evid_Ratio > 39 | Concord_Evid_Ratio < 1/39),
aes(label = gene_symbol),
fill = "white",
size = 3, alpha = 0.9,
min.segment.length = 0, seed = 1,
box.padding = 0.5, point.padding = 0.5,
max.overlaps = 30) +
theme_minimal() +
theme(legend.position = "none",
panel.grid.major.x = element_blank(),
panel.grid.minor.x = element_blank(),
panel.grid.minor.y = element_blank(),
axis.text.x=element_text(size = 15),
#axis.ticks.x=element_blank(),
axis.title.x = element_text(size = 15),
axis.title.y = element_markdown(size = 15))
(Figure_4 <- plot_transcriptome_d/ plot_transcriptome_e + plot_annotation(tag_levels = "a")
)

| Version | Author | Date |
|---|---|---|
| ad196cf | tomkeaney | 2024-04-11 |
Figure 4. Triangles show evidence ratios (ERs) for a \(\Delta R\) (\(R\) on females - \(R\) on males) and b sexual concordance, for all measured traits in the transcriptome dataset. The dashed lines indicate an evidence ratio of 1, where the probability that \(\Delta R\) > 0 is equal to the probability that \(\Delta R\) < 0, or that a trait is just as likely to have a sexually concordant response to selection (ER > 1) as a sexually antagonistic response (ER < 1). The upper dotted line on each plot indicates an evidence ratio of 39, while the lower dotted line indicates an evidence ratio of 1/39; these correspond strongly with the frequentist \(p\) = 0.05. Traits with the most extreme evidence ratios are labelled.
\(~\)
Table S3 / Supplementary dataset
S2: \(R\) partitions and
associated evidence ratios for the gene expression trait dataset. Data
can be downloaded as a .csv file from the github
repository associated with this project.
all_transcriptome_summary <-
R_summarised_transcriptome %>%
select(-c(sd, absolute_R)) %>%
pivot_wider(names_from = c(Fitness_Sex), values_from = 2:4) %>%
rename(`R female` = median_Female,
`R female 2.5% CI` = `2.5%_Female`,
`R female 97.5% CI` = `97.5%_Female`,
`R male` = median_Male,
`R male 2.5% CI` = `2.5%_Male`,
`R male 97.5% CI` = `97.5%_Male`) %>%
left_join(
R_overall_transcript_summarised %>%
select(-sd) %>%
rename( R = R_overall,
`R 2.5% CI` = `2.5%`,
`R 97.5% CI` = `97.5%`),
by = c("Trait", "Trait_Sex")) %>%
left_join(
R_diff_transcript_summarised %>%
select(-sd) %>%
rename(`Delta R` = R_diff,
`Delta R 2.5% CI` = `2.5%`,
`Delta R 97.5% CI` = `97.5%`),
by = c("Trait", "Trait_Sex")) %>%
left_join(evidence_ratios_transcriptome %>%
mutate(Trait_Sex = case_when(Trait_Sex == "Traits measured in females" ~ "Female",
Trait_Sex == "Traits measured in males" ~ "Male")) %>%
select(-c(entrez_id, chromosome, chromosome_numeric, position)) %>%
rename(`R ER` = Overall_Evid_Ratio, `R female ER` = Female_Evid_Ratio,
`R male ER` = Male_Evid_Ratio, `Delta R ER` = Diff_Evid_Ratio,
`Sexual concordance ER` = Concord_Evid_Ratio)) %>%
rename(`Sex trait was measured in` = Trait_Sex) %>%
select(Trait, gene_name, gene_symbol, `Sex trait was measured in`, R, `R 2.5% CI`, `R 97.5% CI`, `R ER`,
`R female`, `R female 2.5% CI`, `R female 97.5% CI`, `R female ER`,
`R male`, `R male 2.5% CI`, `R male 97.5% CI`, `R male ER`,
`Delta R`, `Delta R 2.5% CI`, `Delta R 97.5% CI`, `Delta R ER`, `Sexual concordance ER`) %>%
mutate(across(c(5,6,7,9,10,11,13,14,15,17,18,19), ~round(.x, 3)),
across(ends_with("ER"), ~ round(.x, 4)))
write_csv(all_transcriptome_summary, "data/transcriptome_output/Supplementary_dataset_2.csv")
my_data_table(all_transcriptome_summary)
\(~\)
\(~\)
Fit the model to test whether \(\overline{|R|}\) depends on the sex fitness and trait values were measured in:
# fit the model
median_R_transcriptome_model <-
brm(absolute_R | weights(1/sd) ~ 1 + Fitness_Sex * Trait_Sex + (1|Trait),
family = brmsfamily(family = "Gamma"), # gamma is appropriate for the half-normal distribution created by taking the absolute
data = R_summarised_transcriptome,
prior = c(prior(normal(-2.2, 1), class = Intercept),
prior(exponential(1), class = sd),
prior(exponential(1), class = shape),
prior(normal(0, 0.25), class = b)),
warmup = 2000, iter = 6000,
seed = 1, cores = 4, chains = 4,
control = list(adapt_delta = 0.8, max_treedepth = 10),
file = "fits/median_R_transcriptome_model")
print(median_R_transcriptome_model)
Family: gamma
Links: mu = log; shape = identity
Formula: absolute_R | weights(1/sd) ~ 1 + Fitness_Sex * Trait_Sex + (1 | Trait)
Data: R_summarised_transcriptome (Number of observations: 57072)
Draws: 4 chains, each with iter = 6000; warmup = 2000; thin = 1;
total post-warmup draws = 16000
Group-Level Effects:
~Trait (Number of levels: 14268)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept) 0.41 0.00 0.41 0.42 1.00 2388 5240
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Rhat
Intercept -2.72 0.00 -2.72 -2.71 1.00
Fitness_SexMale -0.01 0.00 -0.01 -0.00 1.00
Trait_SexMale -0.11 0.00 -0.11 -0.10 1.00
Fitness_SexMale:Trait_SexMale 0.03 0.00 0.03 0.04 1.00
Bulk_ESS Tail_ESS
Intercept 1945 4332
Fitness_SexMale 15332 12393
Trait_SexMale 15791 12963
Fitness_SexMale:Trait_SexMale 15058 12576
Family Specific Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
shape 1.65 0.00 1.65 1.66 1.00 23558 11868
Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).
invisible(gc())
To check how does the gamma distribution fares at capturing the shape of the data, we use the model to recapitulate the data it was trained on. Given the shape of the predicted data is similar to the real data, there’s a high chance we’ve used an appropriate distribution family.
pp_check(median_R_transcriptome_model)

| Version | Author | Date |
|---|---|---|
| ad196cf | tomkeaney | 2024-04-11 |
\(~\)
Get model predictions. Once again, we only do this once and save for later convenience.
new_data <- expand_grid(Fitness_Sex = c("Female", "Male"),
Trait_Sex = c("Female", "Male"))
if(!file.exists("data/transcriptome_output/R_transcriptome_fitted.csv")){
R_transcriptome_fitted <-
fitted(median_R_transcriptome_model, newdata = new_data, re_formula = NA, summary = F) %>%
as.data.frame() %>%
rename(FemaleFitness_FemaleTrait = V1, FemaleFitness_MaleTrait = V2,
MaleFitness_FemaleTrait = V3, MaleFitness_MaleTrait = V4) %>%
as_tibble() %>%
mutate(percent_diff_female =
((MaleFitness_FemaleTrait - FemaleFitness_FemaleTrait) / FemaleFitness_FemaleTrait)*100,
percent_diff_male =
((MaleFitness_MaleTrait - FemaleFitness_MaleTrait)/ FemaleFitness_MaleTrait)*100) %>%
pivot_longer(cols = everything(), names_to = "Parameter", values_to = "R_mean") %>%
mutate(Fitness_Sex = case_when(str_detect(Parameter, "FemaleFitness") ~ "Female",
str_detect(Parameter, "MaleFitness") ~ "Male"),
Trait_Sex = case_when(str_detect(Parameter, "FemaleTrait") ~ "Female",
str_detect(Parameter, "MaleTrait") ~ "Male",
str_detect(Parameter, "percent_diff_female") ~ "Female",
str_detect(Parameter, "percent_diff_male") ~ "Male"))
write_csv(R_transcriptome_fitted, "data/transcriptome_output/R_transcriptome_fitted.csv")
} else R_transcriptome_fitted <- read_csv("data/transcriptome_output/R_transcriptome_fitted.csv")
\(~\)
Create panels b and d
R_t1 <-
R_transcriptome_fitted %>%
filter(!str_detect(Parameter, "percent")) %>%
ggplot(aes(x = R_mean, y = Trait_Sex, fill = Fitness_Sex)) + #+ y = Trait_Sex)) + #fill = Fitness_Sex)) +
stat_slab(alpha = 0.9, shape = 21) +#,
labs(x = expression("Mean |"* italic(R) * "| for gene expression traits"),
y = "Phenotyped sex",
fill = "Sex under selection") +
scale_fill_manual(values = c(carto_pal(7, "Purp")[5], carto_pal(7, "Peach")[5])) +
coord_cartesian(xlim = c(0.045, 0.10)) +
theme_minimal() +
theme(panel.background = element_rect(fill='transparent'),
#panel.grid.major.y = element_blank(),
panel.grid.minor.y = element_blank(),
plot.background = element_rect(fill='transparent', color=NA),
legend.position = "bottom",
text = element_text(size=13))
R_t2 <-
R_transcriptome_fitted %>%
filter(str_detect(Parameter, "percent")) %>%
ggplot(aes(x = R_mean, y = Trait_Sex)) + #, #y = Trait_Sex, fill = Trait_Sex)) +
stat_halfeye(.width = c(0.66, 0.95), alpha = 0.9, fill = carto_pal(7, "Emrld")[2],
point_interval = "median_qi", point_fill = "white",
shape = 21, point_size = 4, stroke = 1.2) + # width indicates the uncertainty intervals: here we have 66% and 95% intervals + # width indicates the uncertainty intervals: here we have 66% and 95% intervals+
#scale_fill_manual(values = c(carto_pal(7, "Peach")[5], carto_pal(7, "Peach")[5], carto_pal(7, "Purp")[5])) +
geom_vline(xintercept = 0, linetype = 2, linewidth = 1.2) +
coord_cartesian(xlim = c(-5, 45)) +
xlab(expression("% sex difference in |"*italic(R)* "|")) +
ylab("Phenotyped sex") +
theme_minimal() +
theme(panel.background = element_rect(fill='transparent'),
#panel.grid.major.y = element_blank(),
panel.grid.minor.y = element_blank(),
plot.background = element_rect(fill='transparent', color=NA),
legend.position = "none",
text = element_text(size=13))
Produce the complete figure
left_col <- (R_p1 / R_t1)
right_col <- (R_p2 / R_t2)
wrap_plots(left_col, right_col) &
plot_annotation(tag_levels = "a")

| Version | Author | Date |
|---|---|---|
| ad196cf | tomkeaney | 2024-04-11 |
Figure 2. a the posterior distribution for |\(\overline{R}\)| across organismal level traits, the mean expected response to selection (absolute value) in females and males and split by the sex the trait was measured in. b shows |\(\overline{R}\)| across gene-expression traits. |\(\overline{R}\)| is expressed in trait standard deviations. c and d show the percentage increase or decrease in |\(\overline{R}\)| through males compared to through females for the organismal level and gene expression traits, respectively. Points indicate the estimated mean with associated 66 and 95% credible intervals (where wide enough to be visible).
Calculating stats for the results section
R_transcriptome_fitted %>%
group_by(Parameter) %>%
median_qi(R_mean)
# A tibble: 6 × 7
Parameter R_mean .lower .upper .width .point .interval
<chr> <dbl> <dbl> <dbl> <dbl> <chr> <chr>
1 FemaleFitness_FemaleTrait 0.0661 0.0656 0.0666 0.95 median qi
2 FemaleFitness_MaleTrait 0.0593 0.0588 0.0597 0.95 median qi
3 MaleFitness_FemaleTrait 0.0655 0.0650 0.0660 0.95 median qi
4 MaleFitness_MaleTrait 0.0608 0.0603 0.0612 0.95 median qi
5 percent_diff_female -0.890 -1.47 -0.316 0.95 median qi
6 percent_diff_male 2.55 1.96 3.13 0.95 median qi
sessionInfo()
R version 4.3.1 (2023-06-16 ucrt)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 11 x64 (build 22621)
Matrix products: default
locale:
[1] LC_COLLATE=English_Germany.utf8 LC_CTYPE=English_Germany.utf8
[3] LC_MONETARY=English_Germany.utf8 LC_NUMERIC=C
[5] LC_TIME=English_Germany.utf8
time zone: Europe/Berlin
tzcode source: internal
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] bigsnpr_1.12.4 bigstatsr_1.5.12 ggnewscale_0.4.9 ggrepel_0.9.4
[5] ggtext_0.1.2 broom_1.0.5 bayestestR_0.13.1 tidybayes_3.0.6
[9] brms_2.20.4 Rcpp_1.0.11 patchwork_1.1.3 DT_0.30
[13] kableExtra_1.3.4 rcartocolor_2.1.1 MetBrewer_0.2.0 lubridate_1.9.3
[17] forcats_1.0.0 stringr_1.5.0 dplyr_1.1.3 purrr_1.0.2
[21] readr_2.1.4 tidyr_1.3.0 tibble_3.2.1 ggplot2_3.4.4
[25] tidyverse_2.0.0 workflowr_1.7.1
loaded via a namespace (and not attached):
[1] tensorA_0.36.2 rstudioapi_0.15.0 jsonlite_1.8.7
[4] magrittr_2.0.3 farver_2.1.1 rmarkdown_2.25
[7] fs_1.6.3 vctrs_0.6.4 base64enc_0.1-3
[10] webshot_0.5.5 htmltools_0.5.6.1 distributional_0.3.2
[13] curl_5.1.0 sass_0.4.7 StanHeaders_2.26.28
[16] bslib_0.5.1 htmlwidgets_1.6.2 plyr_1.8.9
[19] zoo_1.8-12 cachem_1.0.8 commonmark_1.9.0
[22] whisker_0.4.1 igraph_1.5.1 iterators_1.0.14
[25] mime_0.12 lifecycle_1.0.4 pkgconfig_2.0.3
[28] colourpicker_1.3.0 Matrix_1.6-1.1 R6_2.5.1
[31] fastmap_1.1.1 shiny_1.7.5.1 digest_0.6.33
[34] colorspace_2.1-0 ps_1.7.5 rprojroot_2.0.3
[37] crosstalk_1.2.0 labeling_0.4.3 fansi_1.0.5
[40] timechange_0.2.0 httr_1.4.7 abind_1.4-5
[43] compiler_4.3.1 rngtools_1.5.2 bit64_4.0.5
[46] doParallel_1.0.17 withr_2.5.1 backports_1.4.1
[49] inline_0.3.19 shinystan_2.6.0 bigassertr_0.1.6
[52] QuickJSR_1.0.7 pkgbuild_1.4.2 gtools_3.9.4
[55] loo_2.6.0 tools_4.3.1 httpuv_1.6.11
[58] threejs_0.3.3 glue_1.6.2 callr_3.7.3
[61] nlme_3.1-162 promises_1.2.1 cmdstanr_0.6.0
[64] bigparallelr_0.3.2 gridtext_0.1.5 grid_4.3.1
[67] checkmate_2.2.0 getPass_0.2-2 reshape2_1.4.4
[70] generics_0.1.3 gtable_0.3.4 bigsparser_0.6.1
[73] tzdb_0.4.0 flock_0.7 data.table_1.14.8
[76] hms_1.1.3 xml2_1.3.5 utf8_1.2.3
[79] foreach_1.5.2 pillar_1.9.0 ggdist_3.3.0
[82] markdown_1.10 vroom_1.6.4 posterior_1.4.1
[85] later_1.3.1 lattice_0.21-8 bit_4.0.5
[88] tidyselect_1.2.0 miniUI_0.1.1.1 knitr_1.44
[91] git2r_0.32.0 arrayhelpers_1.1-0 gridExtra_2.3
[94] V8_4.4.0 svglite_2.1.2 stats4_4.3.1
[97] xfun_0.40 bridgesampling_1.1-2 matrixStats_1.0.0
[100] rstan_2.32.3 stringi_1.7.12 yaml_2.3.7
[103] evaluate_0.22 codetools_0.2-19 cli_3.6.1
[106] RcppParallel_5.1.7 shinythemes_1.2.0 xtable_1.8-4
[109] systemfonts_1.0.5 munsell_0.5.0 processx_3.8.2
[112] jquerylib_0.1.4 coda_0.19-4 svUnit_1.0.6
[115] parallel_4.3.1 rstantools_2.3.1.1 ellipsis_0.3.2
[118] prettyunits_1.2.0 doRNG_1.8.6 dygraphs_1.1.1.6
[121] bayesplot_1.10.0 Brobdingnag_1.2-9 viridisLite_0.4.2
[124] mvtnorm_1.2-3 scales_1.3.0 xts_0.13.1
[127] insight_0.19.6 crayon_1.5.2 rlang_1.1.1
[130] cowplot_1.1.1 rvest_1.0.3 shinyjs_2.1.0